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1.  INTRODUCTION.  The  purpose  of  this  study  was  to  establish  a  baseline 
for  assessing  the  execution  efficiency  of  Model  Based  Vision  (MBV) 
algorithms  coded  in  both  the  Ada  and  C  programming  languages.  Model- 
based  vision  systems  compare  predictions  of  target  signature  in  sensed 
data  with  information  extracted  from  sensed  data  to  achieve  target 
recognition.  To  date,  prediction  of  electromagnetic  sig.iature  for 
Synthetic  Aperture  Radar  (SAR)  is  the  emphasized  sensed  data.  Under 
contract  with  Wright  Laboratory,  The  Analytic  Science  Corporation 
(TASC)  has  extracted  algorithms  from  their  MBV  system  and  delivered  them 
to  the  Government.  The  algorithms  delivered  to  the  Government  are  to  be 
used  as  benchmarks  for  evalutating  the  performance  of  different 
processors  for  MBV  applications.  All  of  the  algorithms  delivered  are 
coded  in  the  C  programming  language. 

To  evaluate  the  suitability  of  Ada  for  MBV  applications,  an  important  MBV 
algorithm  from  the  TASC  MBV  system,  called  Fast  Pairwise  Nearest  Neighbor 
(FASTPNN),  was  translated  in  its  orignal  C  form  to  Ada  and  benchmarked  on 
both  a  VAX  11/780  and  MIPS  Magnum  3000  computer.  While  the  FASTPNN 
algorithm  comprises  only  a  small  subset  of  the  entire  TASC  MBV  system,  in 
terms  of  the  general  types  of  processing  requirements  it  imposes,  it  is 
representative  of  a  significant  portion  of  the  TASC  MBV  system. 


1 


2.  DESCRIPTION  OF  PAIRWISE  NEAREST  NEIGHBOR  ALGORITHM. 


The  Pairwise  Nearest  Neighbor  (PNN)  algorithm  is  a  vector  quantization 
procedure  used  inside  the  information  extraction  portion  of  the  TASC  MBV 
system.  The  PNN  algorithm  is  used  to  represent  the  significant 
characterists  of  a  large  number  of  image  samples  (or  vectors)  with  a 
smaller  specified  number  of  vectors  such  that  the  derived  set  of  vectors 
represent  the  original  set  of  vectors  as  "best"  as  possible  with  respect 
to  coordinate  position  and  weight  (gray  level).  The  PNN  algorithm 
derives  a  specified  number  of  quantization  vectors  by  progressively 
merging  together  pairs  of  vectors  with  minimal  weighted  distance  between 
their  centroids. 

The  TASC  MBV  system  uses  a  vector  quantizer  which  consists  of  the 
FASTPNN  vector  quantizer  cascaded  with  a  Linde-Buzo-Gray  (LBG)  vector 
quantizer.  The  LBG  vector  quantizer  Is  an  iterative  algorithm  that 
produces  a  set  of  quantization  vectors  that  minimize  mean  square 
quantization  error;  however  the  algorithm  tends  to  converge  to  local 
minima  unless  good  estimates  for  initial  quantization  vectors  are 
available.  The  PNN  vector  quantizer  provides  a  near  optimal  set  of 
quantization  vectors,  that  support  the  functional  needs  of  the  LBG 
algorithm. 


2.1  Full  Search  Implementation  of  PNN  Algorithm 

The  key  to  quick  execution  of  the  PNN  algorithm  is  quickly  finding  the 
closest  pairs  of  entries  (vectors)  among  the  distributed  set.  Figure  1 
describes  the  sequence  of  operations  in  the  full  search  implementation 
of  PNN.  Figure  1  shows  that  the  full  search  implementation  of  FNN  is  an 
iterative  procedure  where  on  each  iteration  the  pair  of  vectors  which 
generates  the  least  error  when  merged  (the  closest  pair)  is  merged.  The 
algorithm  explicitly  requires  that,  at  each  iteration,  each  vector's 
nearest  neighbor  be  found.  The  calculation  to  determine  the  closest 
vector  pair  is  essentially  a  weighted  distance  calculation  where  the 
distance  between  the  two  vectors  are  weighted  by  the  gray  level  (or 
weight)  of  each  of  the  two  vectors.  The  new  vector  is  chosen  to 
minimize  the  error  incurred  by  replacing  the  two  original  vectors  into 
one  new  vector. 

In  general,  the  number  of  distance  calculations  to  find  the  two  closest 
vectors  among  a  randomly  distributed  set  is  given  by 

DIST  CALC  =  (N  *  (N-D)  /  2  (1) 

where  DIST_CALC  is  the  number  of  permutations  of  distance  calculations 
required  and  N  is  the  number  of  entries  (vectors).  Thus,  if  one  were 
to  intitiate  a  full  search  PNN  with  1000  entries,  then  the  first 
iteration  would  require  499,500  distance  calculations  to  find  the 
closest  vector  pair  for  merging.  Once  the  first  iteration  was 
completed,  the  following  iterations  would  require  much  fever  distance 
calculations  to  find  the  closest  vector  pairs  since  only  the  distances 
from  the  "new"  vector  to  each  of  the  remaining  entries  would  need  to 
be  recalculated.  The  distance  calculations  among  the  remaining  entries 
would  not  need  to  be  recalculated.  Thus,  in  our  example,  only  999 
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Determine  the  pair  of  vectors  from  data 
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Figure  1 .  Block  Diagram  of  FufI  Search  PNN 


distance  calculations  would  be  required  on  the  second  iteration,  and  998 
iterations  required  on  the  third  iteration.  This  procedure  continues 
Iteratively  until  the  desired  number  of  quantization  vectors  are 
reached.  Thus,  the  total  number  of  distance  calculations  is  given  by 

TOT_^DIST_CALC  -  N*(N-l)/2  +  N-1  +  N-2  +  ...  P  (2) 

where  TOT_DIST_CALC  is  the  total  number  of  distance  calculations 
required,  N  is  the  input  data  size,  and  P  is  the  desired  number  of 
quantization  vectors. 

The  first  term  on  the  right  hand  side  (RHS)  of  equation  (2)  dominates 
the  equation  and  it  is  obvious  that  the  full  search  implementation  of 
PNN  has  computational  requirement  proportional  to  0(N**2). 


2.2  Fast  Implementation  of  PNN  Algorithm 

If  the  requirement  to  merge  the  absolute  closest  pair  of  vectors  at 
each  step  is  relaxed,  as  long  as  vectors  get  merged  eventually  then  a 
fast  implementation  of  the  PNN  algorithm  can  be  derived  (Reference  2). 
This  approximated  PNN  algorithm  is  called  the  "fast"  PNN  (FASTPNN) 
algorithm.  For  most  applications  the  error  introduced  by  the  FASTPNN 
approximation  of  PNN  is  considered  tolerable  (reference).  The  "fast" 
Implementation  of  PNN  was  used  in  this  study. 

Figure  2  describes  the  sequence  of  operations  in  FASTPNN.  The  FASTPNN 
algorithm  initially  constructs  a  set  of  entries  from  the  input  data, 
then  iteratively  merges  pairs  of  entries  until  a  specified  number  of 
quantization  vectors  are  reached.  At  each  iteration,  the  algorithm 
proceeds  by  recursively  splitting  sets  of  entries,  or  buckets,  into 
pairs  of  buckets  until  the  total  number  of  entries  in  each  bucket  is 
at  or  below  a  prespecified  number.  From  each  bucket,  a  candidate  pair 
of  entries  is  nominated  that,  if  merged,  provide  the  smallest  increase 
in  quantization  error  of  all  pairs  within  the  bucket.  A  specified 
percentage  of  the  pairs  are  selected  and  merged,  with  the  individual 
buckets  merged  into  a  single  bucket  to  complete  an  iteration.  The 
sequence  is  repeated  until  the  total  number  of  entries  is  equal  to  the 
specified  number  of  quantization  vectors.  The  value  vectors  in  the 
final  set  of  entries  define  the  derived  quantization  vectors. 

The  computational  savings  in  FASTPNN  is  derived  in  the  algorithm's 
ability  to  "intelligently"  split  the  global  data  set  among  a  specified 
number  of  localized  regions  or  buckets.  Merging  of  vectors  are  then 
performed  independently  and  with  more  efficiency  within  these  localized 
regions  (buckets). 

A  key  parameter  in  FASTPNN  is  BUCKETSIZE,  the  maximum  number  of  vectors 
per  bucket.  BUCKETSIZE  plays  a  key  role  in  determining  the  number  of 
distance  calculations  involved  in  finding  the  closest  vector  pair. 

Since  FASTPNN  requires  that  the  data  be  split  evenly  across  all  buckets 
then  the  total  number  of  buckets  required  for  a  particular  iteration 
equals  the  number  of  input  vectors  divided  by  BUCKETSIZE. 

BUCKET  NUM  -  N/BUCKETSIZE  (3) 
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Block  1  Organize  data  samples  (veclofs)  into  a  K-d  tree 
containing  one  bucket _ _ 
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Figure  2.  Block  Diagram  of  FASTPNN  Algorithm 


where  BUCKET_NUM  equals  the  total  number  of  buckets  for  a  data  set  of 
size  N,  and  bucket  containing  BUCKETSIZE  entries.  From  equation  (1),  to 
find  the  closest  vector  pair  in  each  particular  bucket  requires  that 
(BUCKETSIZE  *  (BUCKETSIZE  -  1))  /  2  weighted  distance  calculations  be 
performed.  Thus,  the  total  number  of  distance  calculations  for  the 
complete  data  set  on  tlie  first  iteration  equals  the  number  of  distance 
calculations  per  bucket  times  the  number  of  buckets,  or 

DISTCALC  =  (BUCKETSIZE  *  (BUCKETSIZE  -  l)/2)  *  N/BUCKETSIZE  (4) 

Equation  (4)  reduces  to 

DIST_CALC  =  N  *  (BUCKETSIZE  -  l)/2  (5) 

where  DIST_CALC  is  the  total  number  of  distance  calculations  required  on 
the  first  iteration  of  FASTPNN.  Thus,  it  is  readily  seen  from  equation 
5  that  as  BUCKETSIZE  decreases  the  number  of  calculations  goes 
proportionally  down.  In  the  "worst-case"  limiting  case,  if  BUCKETSIZE 
is  equal  to  N  (all  of  the  entries  are  in  one  bucket)  we  get  the  same 
result  as  obtained  in  equation  (1)  for  the  full  search  implementation. 

As  a  simplified  example  of  the  magnitude  by  which  FASTPNN  reduces  the 
number  of  distance  calculations,  let's  suppose  as  we  did  in  the  last 
section  that  we  begin  PNN  with  1000  entries.  If  we  assume  that 
BUCKETSIZE  is  ten,  then  only  4,500  distance  calculations  are  required  on 
the  first  iteration  to  find  the  closest  vector  pairs.  This  is 
significantly  lower  than  the  499,500  distance  calculations  required  by 
the  full  search  implementation  on  the  first  iteration. 

The  number  of  vectors  merged  per  iteration  is  given  by  the  term 
KDMERGE  *  BUCKETNUM  where  KDMERGE  is  a  fixed  percentate  of  the  top 
vector  pair  candidate  from  each  bucket. 

Thus,  on  the  second  iteration,  the  number  of  distance  calculations  will 
be  reduced  to  (N  -  KDMERGE*N/BUCKF.TSIZE)  *  (BUCKETSIZE  -1)  /2 

Although,  the  number  of  distance  calculations  in  determining  the 
candidate  vector  pairs  for  merging  are  greatly  reduced  by  using  FASTPNN, 
there  are  other  computational  "overhead"  factors  associated  with  the 
FASTPNN  algorithm  which  are  not  present  in  the  full  search 
implementation  of  PNN.  From  Figure  2,  these  include  Initialization  of 
the  K-d  tree  data  structure  and  splitting  and  merging  of  buckets. 
Despite  these  additional  "overhead"  factors,  it  has  been  shown  that 
FASTPNN  has  computational  efficiency  proportional  to  0(N  log  N). 


6 


3.  SOFTWARE/BENCHMARK  IMPLEMENTATION  OF  FASTPNN 


As  stated  in  the  introduction,  the  FASTPNN  algorithm  was  delivered  to 
the  Government  coded  in  the  C  programming  language.  Figure  3  shows  a 
subprogram  calling  tree  of  the  FASTPNN  algorithm.  Note  from  Figure  3, 
that  the  FASTPNN  subprogram  structure  divides  into  3  main  parts.  These 
three  parts  include  BuildRDtree,  MergedownKDtree,  and  DestroyKDtree. 

BuildKDtree  contains  the  subroutines  which  involve  building  and 
initialization  of  the  K-d  tree  data  structure.  The  K-d  tree  is  the 
data  structure  which  permits  the  partitioning  of  the  global  data  set 
into  buckets  (localized  regions)  from  which  nearest  neighbor  searches 
can  be  performed  independently.  With  reference  to  Figure  2,  the 
functionallity  of  BuildKDtree  is  described  by  block  1  of  of  the  block 
diagram. 

MergeDownKDtree  contains  the  routines  which  reduce  the  K-d  tree  built  in 
Buildkdtree  to  the  specified  number  of  centroids  (vectors).  With 
reference  to  Figure  2,  the  functionallity  of  MergeDownKDtree  is 
described  by  blocks  two  through  six  of  the  block  diagram. 

DestroyKDtree  contains  the  routines  which  permit  the  deallocation  of  memory 
which  was  dynamically  allocated  in  BuildKDtree  and  MergeDownKDtree.  The 
functionallity  of  DestroyKDtree  is  not  described  in  Figure  2,  as  it  is 
not  a  necessary  part  of  the  FASTPNN  algorithm.  DestroyKDtree  is  merely 
provided  to  implement  the  good  software  practice  of  deallocating  objects 
no  longer  in  use. 


3.1  C  Implementation  Structure  of  FASTPNN  Benchmark 

To  permit  timing  of  the  C  coded  FASTPNN  algorithm,  FASTPNN  was 
integrated  with  a  timer  to  form  a  benchmark.  The  timer  measured  the 
process  CPU  time  utilized  while  running  FASTPNN.  The  basic  strategy  in 
timing  the  FASTPNN  algorithm  was  to  have  a  main  (driver)  program  obtain 
the  initial  time  just  before  making  a  function  call  to  FASTPNN,  followed 
by  another  call  to  the  timer  just  prior  to  executing  FASTPNN.  The 
difference  between  these  two  times  is  the  time  of  interest  (elapsed 
time).  Because  timers  use  system  dependent  resources,  a  separate  timer 
function  was  written  for  both  the  VAX  11/780  and  the  MIPS  MAGNUM  3000. 
Appendix  A  provides  the  C  source  code  listing  of  the  FASTPNN  benchmark. 

Figure  4  shows  the  module  structure  of  the  C  coded  FASTPNN  benchmark. 

The  FASTPNN  benchmark  is  composed  of  five  modules.  Two  of  the  five 
modules  are  header  files.  Header  files  are  used  in  C  to  contain 
definitions  and  declarations  which  are  to  be  shared  among  different 
files  (modules).  The  arrows  in  Figure  4  are  used  to  indicate  the 
dependency  of  particular  modules  on  the  header  files.  For  instance,  the 
arrow  from  FASTPNN. c  to  FASTPNN. h  indicates  that  FASTPNN. h  is  to  be 
"included"  into  FASTPNN. c.  Header  file  Timer. h  contains  the  interface  to 
the  system  routines  which  perform  the  timing  of  the  benchmark.  Header 
file  FASTPNN. h  contains  the  data  structure  declarations  which  are  used 
in  FASTPNN. c  and,  to  a  limited  extent,  in  Main.c.  The  main  program, 
Main.c,  is  used  to  read  the  data,  perform  the  timing  of  FASTPNN,  and 
output  the  results. 
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Figure  3.  FASTPNN  Subprogram  Calling  Tree 


FASTPNN.c 


3.2  Ada  Implementation  of  FASTPNN  Benchmark 


Presented  in  this  section  is  the  organization  structure  from  a 
module  (or  compilation  unit)  level  of  the  C  to  Ada  translated  FASTPNN 
benchmark.  The  details  of  the  actual  coding  from  C  to  Ada  are  addressed 
in  section  4  of  this  report.  Appendix  B  provides  the  source  code 
listing  of  the  Ada  coded  FASTPNN  benchmark. 

The  Ada  coded  FASTPNN  algorithm  was  integrated  with  a  timer  to  form 
a  benchmark.  The  timer  measured  the  process  CPU  time  utilized  while 
running  FASTPNN.  As  with  the  case  of  the  C  benchmark,  a  separate  timer 
(body)  was  written  for  both  the  VAX  11/780  and  the  MIPS  MAGNUM  3000. 

Figure  5  shows  the  overall  module  structure  of  the  Ada  coded  FASTPNN 
benchmark.  Five  packages  were  used  in  the  benchmark  implementation. 
Data_Struct_Pkg  contains  the  declarations  of  all  data  types  used  in  the 
benchmark.  Packages  Build_Pkg,  MergeDownPkg ,  and  Destroy_Pkg  contain 
all  of  the  FASTPNN  subroutines  presented  in  Figure  3.  Package  Timer_Pkg 
contains  the  sytem  dependent  routines  which  permit  benchmark  timing  to 
be  performed.  The  arrows  in  Figure  5  are  used  to  indicate  the  dependency 
of  the  particular  modules  with  each  other.  For  instance,  the  arrows 
pointing  from  FASTPNN  to  Build  Pkg,  MergeDown_Pkg,  and  Destroy_Pkg 
indicates  that  these  three  packages  must  be  "withed  into"  Build_Pkg. 

The  main  program  is  used  to  input  the  data,  time  FASTPNN,  and  output  the 
results. 

Figure  6  shows  the  "package  location"  of  each  of  the  FASTPNN 
subroutines.  Build_Pkg  contains  all  of  the  routines  which  are  nested  in 
the  Bui Idkd tree  portion  of  Figure  3.  MergeDown_Pkg  contains  the 
routines  which  are  nested  in  the  MergeDownKDtree  portion  of  Figure  3. 
DestroyPkg  contains  the  routines  which  are  nested  in  the  DestroyKdtree 
portion  of  Figure  3.  Listed  in  the  Ada  specification  portion  of  each  of 
the  above  three  packages,  are  the  routines  which  must  be  visible  outside 
of  their  local  package  usage.  Note  from  Figure  3,  that  since  FASTPNN 
calls  procedures  BuildKDtree,  MergeDownKDtree,  and  DestroyKDtree,  these 
three  procedures  are  all  included  in  the  specification  (as  opposed  to 
body)  portion  of  their  respective  packages.  Procedure  CreateKDbucket  is 
included  in  the  specification  portion  of  Build_Pkg,  since  it  called  by  a 
routine  (SplitBucket)  outside  of  Build_Pkg.  Procedures  DestroyKDnode, 
DestroyKDbucket,  and  DestroyKDentry  are  included  in  the  specification 
portion  of  Destroy_Pkg  since  they  are  called  by  routines  outside  of 
DestroyPkg. 


3.3  Benchmark  Data  Set 

TASC  provided  the  Government  with  an  input  data  set  to  be  used  for 
executing  the  FASTPNN  benchmark.  The  input  data  set  consisted  of  a 
single  chip  of  SAR  signature  prediction  data  of  a  B52  aircraft.  A  chip 
is  a  subregion  of  interst  extracted  from  a  larger  image  (i.e  IK  x  IK). 
Chips  typically  range  in  size  fro.i  64x64  to  128x128  pixels.  An  image 
plot  of  the  B52  aircraft  chip  is  displayed  in  Figure  7.  This  data  was 
obtained  by  using  a  prediction  tool  developed  at  TASC  for  predicting  SAR 
returns  from  various  targets. 
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Figure  5.  Ada  Module  Breakdown  of  FASTPNN  Benchmai 


SPECIFICATION 


Figure  6.  Ada  Package  Implementation  of  FASTPNN 


The  input  SAR  data  set  (chip)  contained  96x96  samples  (vectors).  Each 
sample  contained  3  components*  an  x  coordinate  position,  a  y  coordinate 
position,  and  a  gray  level  (weight).  Thresholding  was  used  to  reduce  the 
intial  data  set  o£  5184  vectors  to  1,675  vectors.  The  purpose  of 
thresholding  is  to  reduce  the  data  set  by  eliminating  vectors  with  non 
significant  weight.  The  reduced  set  of  1,675  vectors  was  the  actual 
input  to  the  FASTPNN  algorithm. 
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4.  ADA  TO  C  TRANSLATION  STRATEGY. 

To  establish  a  uniform  framework  tor  recoding  the  FASTPNN  benchmark 
from  C  to  Ada,  a  C  to  Ada  translation  strategy  was  adopted.  Because  the 
goal  is  to  compare  the  inherent  efficiency  of  Ada  with  C,  no  attempt  was 
made  in  the  translation  process  to  improve  the  execution  efficiency  of 
the  Ada  source  code  beyond  that  of  the  C  source  code.  Thus,  a  direct 
"line  by  line"  translation  strategy  was  adopted. 

4.1  Data  Structures 

The  first  step  in  the  translation  process  was  to  "map"  the  individual  C 
data  structures  to  "equivalent"  Ada  data  structures.  In  general,  the 
mapping  between  C  and  Ada  data  stuctures  is  relatively  straightforward. 
Table  1  below  highlights  the  main  data  structure  mappings  used  in  this 
study. 


TABLE  1.  MAPPING  OF  C  TO  ADA  DATA  STRUCTURES 


C 


Ada 


Dynamic  Array  (*) 

Pointer  Type 

Dynamic  Array  of  Pointers  (**) 

Struct 

Union 


Unconstrained  Array 
Access  Type 

Unconstrained  Array  of  Access  Types 
Record 

Record  with  Variant  Part 


Table  2  provides  an  example  of  the  data  structure  mapping  techniques 
presented  in  Table  1  by  presenting  two  of  the  translated  data  structures 
which  were  extracted  from  the  FASTPNN  algorithm.  The  left  hand  side  of  Table 
2  shows  the  original  C  coded  data  structure  and  the  right  hand  side  of  Table 
2  shows  the  Ada  translated  data  stucture.  These  data  structures  can  be  found 
in  Fa.stpnn.h  (Appendix  A)  and  Data_Struct_Pkg  (Appendix  B)  of  the  C  and  Ada 
source  code  respectively. 


The  sections  below  provide  a  brief  explanation  of  the  data  structure 
mappings  displayed  i.i  Table  1  and  used  as  examples  in  Table  2. 


4.1.1  Dynamic  Arrays/Unconstrained  Arrays 

C  Implementation.  Structure  componets  *mean,  *wmean,  and  *wsqmn  of 
struct  kdentry  in  Table  2a  all  correspond  to  dynamic  arrays  containing 
floating  point  numbers.  The  number  of  array  components  in  each  of  the 
above  arrays  is  constrained  at  run  time  to  equal  the  dimension  of  the 
problem's  application  (i.e.  if  the  problem  dimension  is  two,  there  would 
he  two  components  for  each  array!.  Thus,  in  the  case  of  FASTPNN,  the 
use  of  dynamic  arrays  permits  the  fiexiblity  to  use  FASTPNN  for  a 
generalized  n  dimensional  application.  At  run-time,  memory  space  is 
dynamically  allocated  by  making  a  call  to  the  libray  function  calloc 
with  the  dimension  size  passed  as  a  parameter  to  the  function  calloc. 

Ada  Implementation.  In  Ada,  one  also  could  also  dynamically  allocate 
the  arrays  *mean,  ^weight,  and  *wsqmn  by  using  allocators.  However,  Ada 
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TABLE  2.  EXAMPLE  OF  C  TO  ADA  DATA  STRUCTURE 
MAPPING 

Table  2  a 


•nd  ea««: 

•nd  raeord; 

norb-R  :  kdalaaifUnoda-lypa): 
norb  :  hdalbnfkdbuelibt.typb); 


permits  an  easier  and  more  efficient  method  to  derive  most  of  the 
flexibility  of  the  dynamic  array  through  use  of  the  unconstrained  array. 
The  unconstrained  array  enables  one  to  treat  arrays  which 
have  the  same  characteristics  but  differ  only  in  their  sis^e,  as 
equivalent  types.  Type  mean_array_type  in  Table  2a  is  an  example  of  an 
unconstrained  array  type. 

4.1.2  Pointer  Types/Access  Types 

Pointer  types  in  C  map  to  access  types  in  Ada. 


4.1.3  Dynamic  Array  of  Pointers/Unconstrained  Array  of  Access  Types 

C  Implementation.  Dynamic  arrays  of  pointers  were  used  to  maintain  linked 
list  pointers  in  multiple  dimensions.  As  explained  in  section  4.1.1, 
the  size  of  the  arrays  are  dynamically  allocated  to  equal  the  dimension 
of  the  application  problem.  In  Table  2a,  structure  component  **next 
provides  an  example  of  a  dynamic  array  of  pointers. 

Ada  Implementation.  Unconstrained  arrays  of  access  types  were  used  to 
maintain  linked  list  pointers  in  multiple  directions.  In  Table  2a, 
record  componet  next_array  provides  an  example  of  an  unconstrained  array 
of  access  types. 

4.1.4  Structures/Records 

Structures  in  C  map  directly  to  records  in  Ada. 


4.1.5  Structures  Containing  Unions/Variant  Records. 

C  Implementation.  Unions,  like  structures,  contain  members  whose 
individual  data  types  may  differ  from  one  another.  However  the  members 
that  compose  a  union  all  share  the  same  storage  area  within  the 
computer's  memory,  whereas  each  member  within  a  structure  is  assigned 
its  own  unique  storage  area.  Table  2b  provides  an  example  of  a  structure 
containing  a  union. 

ADA  Implementation.  Ada  does  not  contain  an  equivalent  data  structure 
directly  mapable  to  the  C  union  isolated  (not  contained  inside  another 
data  structure).  However,  the  Ada  variant  record  maps  closely  to  the 
case  where  a  where  a  C  union  is  contained  within  a  C  structure.  In  this 
case,  the  "union"  portion  is  represented  by  the  variant  part  of  the 
record.  Table  2b  provides  an  example  of  a  variant  record. 
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^.2  C  to  Ada  Coding  Dualities 


The  final  step  in  the  C  to  Ada  coding  strategy  was  to  highlight  the 
significant  coding  dualities  between  C  and  Ada.  Table  3  below  outlines 
a  non  exhaustive  list  of  coding  dualities  between  C  and  Ada. 

TABLE  3.  Miscellaneous  C  to  Ada  Issues 

C  Ada 

Loop  Statements  Loop  Statements 

-  for  -  for 

-  while  -  while 


Break 


Exit 


Conditional  Statements 

-  if 

-  if/else 

-  if/else  if/else  if/... /else 

-  switch 
Functions 


Conditional  Statements 

-  if/end  if 

-  if/else/end  if 

-  if /elsi f/elsif/ ... /else/end  if 

-  case 

Procedures  or  Functions 


Pointer  Type  Operations 

-  *P 

-  p->struciure  component 

-  &p 


Access  Types 

-  p.all 

-  p. record  componet 

-  p' ADDRESS 


Dynamic  Allocation 
-  calloc,  malloc 


Dynamic  Allocation 
-  new 


Deallocation 
-  Cfree 


Deallocation 

-  Unchecked  Deallocation 


Logical  Operators 
-  && 


Logical  Operators 

-  and  then 

-  or  else 


4.2.1  Loop  Statements 

Loop  statements  map  very  closely  between  C  and  Ada  and  only  Involve 
straightforward  syntax  modifications  to  convert  from  one  language  to  the 
other 

4.2.2  Conditional  Statements 

Conditional  statements  map  very  closely  between  C  and  Ada  and  only 
involve  straightforward  syntax  modifications  to  convert  from  one 
language  to  the  other 

4.2.3  "Break'V'Exi t"  Statement 
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The  "break"  statement  in  C  and  the  "exit"  statement  in  Ada  permit  the 
innermost  enclosing  loop  to  be  exited  immediately  without  "testing"  at 
the  bottom  or  top  of  the  loop. 

4. 2. A  Subprograms 

C  Implementation.  All  subprograms  in  C  are  implemented  as  functions. 

A  C  function  returns  a  single  value  as  a  result  of  a  call. 

Ada  Implementation.  A  subprogram  in  Ada  can  be  implemented  either  as  a 
function  or  a  procedure.  The  following  guidelines  were  used  to  determine 
whether  a  function  or  a  procedui.c  was  used.  A  procedure  was  used  whenever 
two  or  more  parameters  were  to  be  modified  through  a  subroutine  call.  If 
one  or  less  parameters  was  to  be  modified  then  either  a  function  or  a 
procedure  was  used.  The  decision  to  choose  either  a  function  or  a 
procedure  was  based  upon  the  characteristics  of  the  subprogram.  The 
following  coding  convention  was  used  when  deciding  upon  a  function  or  a 
procedure.  It  the  state  of  the  input  argument  was  to  be  modified  then  a 
procedure  call  was  made.  Thus,  if  the  subprogram  could  be  written  as  an 
Ada  "in-out"  variable,  then  an  Ada  procedure  was  used.  If  just  some 
value  was  returned,  then  a  function  was  used.  Thus,  if  the  variable  that 
needed  to  change,  was  characteristic  of  an  Ada  "out"  parameter,  then  an 
Ada  function  was  used. 

4.2.5  Pointer  Operations 

C  Implementation.  The  unary  operator  *  is  called  the  indirection 
operator.  When  applied  to  a  pointer,  the  indirection  operator  accesses 
the  object  that  the  pointer  points  to.  For  example,  if  a  pointer- 
variable,  p,  points  to  an  integer  of  value  2,  then  the  statement  x  ■  *p 
assigns  a  value  of  2  to  the  integer  variable  x.  If  we  assume  that  the 
pointer  variable,  p,  points  to  a  structure  which  contains  two  members 
(structure  components),  say  member_l  and  meraber_2  respectively,  then  the 
statement  x  =  p  ->  member_l  assigns  the  value  of  structure  component 
memberl  to  the  variable  x,  and  the  statement  y  =  p  ->  member_2  assigns 
the  value  of  structure  component  member  2  to  variable  y.  The  unary 
operator  when  applied  to  a  variable,  gives  the  memory  address  of  the 
object.  For  example,  the  statement  p  s  fiix  assigns  the  address  of  x  to 
the  pointer  variable  p. 

Ada  Implementation.  The  word  "all"  is  used  to  access  objects  pointed  to 
by  an  access  type  variable.  For  example,  if  an  access  variable  type,  p, 
points  to  an  integer  of  value  2,  then  the  statement  x  :=  p.all  assigns  a 
value  of  2  to  the  integer  variable  x.  If  we  assume  that  the  access 
variable,  p,  points  to  a  record  which  contains  two  members,  say  membei_l 
and  member_2  respectively,  then  the  statement  x  :=  p.memberl  assigns 
the  value  of  record  component  memberl  to  the  variable  x  and  the 
statement  y  :=  p.member2  assigns  the  value  of  record  componet  member_2 
to  the  variable  y.  To  obtain  the  memory  address  of  a  variable,  the 
attribute  'Address  from  package  system  must  be  used. 

4.2.6  Dynamic  Allocation 

C  Implementation.  The  function  calloc  was  used  in  conjunction  with  the 
function  sizeof  to  obtain  blocks  of  memory  dynamically.  As  an  example, 
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the  statement 


ip  =  calloc(n,  sizGof(int)) 

allocates  memory  for  an  array  having  n  elements,  vith  each  component 
having  a  length  equal  to  number  of  bytes  used  by  the  host  machine  in 
representing  an  integer  (usually  2).  A  pointer  ip  is  returned  to  these 
n*2  (assuming  an  integer  is  2  bytes)  bytes  of  uninitialized  storage,  or 
NULL  if  the  request  cannot  be  satisfied. 

Ada  Implementation.  Dynamic  allocation  is  performed  in  conjunction  with 
the  "new"  statement.  Below  is  the  code  to  translate  the  same  example  as 
presented  in  the  C  case  above  into  Ada 

type  integei_array_type  is  array  (1  ..  n)  of  integer 
ip  ;=  new  intGger_array_type; 

The  first  statement  of  the  above  code  declares  an  array  of  integers  of  n 
elements  long  of  type  integer_array_type.  The  second  statement  uses  the 
keyword  "new"  to  create  a  designated  object  of  type  integer_array_type 
and  it  stores  the  memory  address  of  the  newly  created  object  in  access 
type  ip. 

4.2.7  Memory  Deallocation 

Deallocation  is  used  to  "free"  memory  which  was  previously  dynamically 
allocated  but  is  no  longer  in  use.  The  storage  used  by  that  variable 
can  then  be  reused  to  allocate  a  new  variable. 

C  Implementation.  The  function  free  deallocates  space  which  was 
previously  dynamically  allocated  by  a  call  to  calloc  or  malloc.  For 
example,  if  we  assume  that  memory  was  previously  dynamically 
allocated  by  the  statement  p  =  (int  *)  calloc(n,  sizeof(int))  hen  the 
statement  free(p)  frees  the  space  pointed  to  by  p. 

ADA  Implementation.  UncheckedDeallocation  is  one  of  the  four 
predefined  generic  units  that  are  provided  by  every  implementation  of 
the  Ada  language  .  In  the  Ada  language,  allocated  variables  are 
deallocated  by  calling  an  instance  of  the  generic  procedure 
Unchecked_Deallocation.  A  call  on  an  instance  of  Unchecked_Deallocation 
takes  one  variable,  which  is  a  variable  of  access  type.  The  allocated 
variable  designated  by  the  parameter  is  deallocated,  and  the  parameter  is 
set  to  NULL.  For  example,  assume  that  memory  was  previously  dynamically 
allocated  by  the  statement  p  ;=  new  Integer_Array_Type;  where 
IntGger_Array_Type  is  an  array  of  integers.  To  deallocate  the  array  of 
integers,  it  is  first  necessary  to  instantiate  the  generic  procedure 
UncheckedDeallocation  by  the  statement  "procedure  FREE  is  new 
Unchecked_Deallocation  (Integer,  Integer_Array_Ptr_Type) ; "  and  then  the 
statement  "FREE(Integer_Array_Ptr)"  is  used  deallocate  the  array.  Note, 
that  there  is  no  requirement  to  use  the  name  "FREE"  as  any  name  (i.*^ 
DeallocateCell ,  Dispose,  etc)  could  be  used. 

4.2.8  Logical  Operators. 

C  Implementation.  Logical  operators  joined  by  the  logical  operators 
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&&  and  1  I  are  evaluated  lef t-to-right,  but  only  until  the  overall 
true/false  has  been  established.  In  the  case  of  the  &&  operator,  if  the 
left-hand  expression  evaluates  to  false,  the  AND  THEN  operator  returns 
the  value  false  immediately  without  evaluating  the  right-hand 
expression.  In  the  case  of  the  ] [  operator,  if  the  left-hand  expression 
is  true,  the  | |  operator  returns  a  true  value  without  checking  the  other 
operand . 

ADA  Implementation.  The  Ada  equivalents  of  &&  and  1 |  respectively 
are  AND  THEN  and  OR  ELSE  repectively.  These  are  called  short-circuit 
AND  and  short-circuit  OR  respectively.  They  produce  the  same  results  as 
the  plain  AND  and  OR,  but  they  force  the  computer  to  evaluate  the 
expression  in  lef t-to-right  order. 
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5.  BENCHMARK  DEVELOPMENT  SYSTEMS 

This  section  provides  a  description  of  the  two  development  systems  used 
in  this  study  to  perform  the  benchmarking. 

5.1  MIPS  MAGNUM  3000. 

5.1.1  MIPS  MAGNUM  3000  Development  System.  Below  is  a  description  in 
outline  form  of  the  MIPS  Magnum  3000  workstation  used  in  this 
benchmarking  study 

-  R3000  Central  Processing  Unit  (CPU)  running  at  25  MHZ 

-  R3010  Floating-Point  Coprocessor  (FPC),  running  at  25  MHZ 

32  Kbytes  of  instruction  cache  (I-cache)  and  32  Kbytes  of  data 
cache  (D-cache),  using  20  ns  static  RAM  (SRAM)  with  fixed  8-word 
block-refill  size 

-  ASIC  Read/Write  Buffer  with  8-vord  buffering 

16  Mbytes  of  100  ns  DRAM  which  supports  block-mode  transfers  with 
peak  data  rates  of  50  Mbytes  per  second  on  writes  and  100  Mbytes 
per  second  on  reads 

The  R3000  CPU  provides  32  general  purpose  32-bit  registers,  a  32-bit 
Program  Counter,  and  two  32-bit  registers  that  hold  the  results  of 
integer  divide  and  multiply  operations.  The  R3010  FPC  is  tightly 
coupled  to  the  R3000  CPU  and  can  execute  instructions  in  parallel  with 
the  CPU.  The  R3010  contains  sixteen,  64-bit  registers  that  can  be  used 
to  hold  single-precision  or  double-precision  values.  The  MIPS  R3000  has 
74  instructions  while  the  MIPS  R3010  has  20  instructions.  However, 
since  eight  of  the  R3000  instructions  are  common  to  R3010  instructions, 
there  are  a  combined  total  of  86  MIPS  instructions.  The  R3000  can 
access  memory  only  through  simple  load/store  operations.  All  MIPS 
instructions  are  32  bits  long.  Although  there  is  only  a  single 
addressing  mode  (base  register  plus  16-bit  signed  displacement),  there 
are  numerous  individual  load  and  store  instructions  that  can  load  or 
store  integer  data  in  sizes  8,  16,  and  32  bits  with  signed  and  unsigned 
extension. 

5.1.2  MIPS  Compilers.  MIPS  compilers  support  six  programming  languages 
including  C  and  Ada.  The  compiler  system  has  a  separate  front-end  to 
translate  each  language  and  a  common  back-end  to  generate  optimized 
machine  code.  Run-time  libraries  provide  language-dependent  functions 
for  each  language.  The  front-ends  translate  the  semantics  of  each 
language  into  an  intermediate  representation,  called  U-code.  U-code  is 
used  by  several  of  the  common  back-end  components. 

MIPS  uses  a  common  global  optimizer,  called  uopt,  for  all  of  their 
compilers.  MIPS  categorizes  there  compiler  optimizations  into  four 
different  levels.  Below  is  a  summary  of  the  different  optimization 
levels : 

Level  1  -  includes  peephole  and  local  optimizations 
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Level  2  -  includes  level  one  optimizations,  plus  the  following  global 
optimizations:  loop-invariant  code  motion,  strength  reduction,  common 
subexpression  elimination,  and  register  allocation 

Level  3  -  includes  level  two  optimizations  plus  interprocedure 
register  allocation 

Level  4  -  includes  level  3  optimizations,  plus  procedure  inlining 


5. 1.2.1  MIPS  Ada  Compiler. 

The  MIPS  Ada  compiler  provides  two  global  optimizers,  namely  uopt 
(discussed  in  section  5.1.2)  and  0PTIM3. 

0PTIM3  is  a  high  le/el  global  optimizer  that  performs  many  classical 
code  optimizations  and  several  that  are  specific  to  Ada.  These  include 
redundant  range  check  elimination  and  range  propagation  for  elimination 
of  constraint  checking. 

Version  3.0  of  the  MIPS  Ada  compiler  was  used  in  this  benchmarking 
study.  MIPS  Ada  3,0  is  intended  to  support  the  level  2  optimizations 
presented  in  section  5.1.2.  However,  although  the  Ada  code  compiled 
properly  using  level  2  optimizations,  the  run-time  execution  terminated 
due  to  a  "segmentation  error".  Thus,  it  was  necessary  to  only  use  level 
1  optimizations  in  compiling  the  Ada  source  code. 


5, 1.2. 2  MIPS  C  Compiler. 

The  MIPS  C  compiler  used  was  the  2.11  release.  The  2.11  release  supports 
all  of  the  level  4  optimizations  presented  in  section  5.1.2.  In  this 
study,  the  MIPS  C  code  was  compiled  and  executed  using  level  4 
optimizations . 


5.2  VAX  11/780. 

5.2,1  VAX  Hardware. 

Below  is  a  description  in  outline  forta  which  summarizes  the  main 
features  of  the  VAX  11/780  computer. 

contains  a  VAX  11/780  processor 

-  Floating  Point  Accelerator  is  optional  (used  in  this  study) 

contains  an  8  Kbyte  cache  which  hold  both  data  and  instructions 

contains  as  addresss  translation  buffer  (cache)  which  can  hold  up  to 
128  virtual-to-physical  page-address  translations 

-  contains  an  8  byte  instruction  buffer 
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The  VAX  is  the  classical  example  of  a  Complex  Instruction  Set  Computer 
(CISC)  architecture.  There  are  over  200  different  instructions  and  7 
basic  addressing  modes.  The  instruction  set  operates  on  integer, 
floating-point,  character-strings  and  packed-decimal  strings,  and  bit 
fields.  The  processor  provides  64-bit,  32-bit,  16-bit,  and  8-bit 
arithmetic;  instruction  prefetch;  and  an  address  translation  buffer. 

The  CPU  includes  16  32-bit  general  purpose  registers  for  data 
manipulatuion  and  the  Processor  Status  Longword  for  controlling  the 
execution  states  of  the  CPU.  The  VAX  used  in  this  study  contained  an 
optional  high  performance  floating  point  accelerator  (FPA).  The  FPA  is 
an  independent  processor  that  executes  in  parallel  with  the  base  CPU. 
The  FPA  takes  advantage  of  the  CPU's  instruction  buffer  to  access  main 
memory.  Once  the  CPU  has  the  required  data,  the  FPA  overrides  the 
normal  execution  flow  of  the  standard  floating-point  microcode  and 
forces  use  of  its  own  code.  While  the  FPA  is  executing  the  CPU  can  be 
performing  other  operations  in  parallel. 


5.2.2.  VAX  COMPILERS.  This  section  describes  both  the  VAX  Ada  compile  and 
also  the  VAX  C  compiler, 

5. 2. 2.1  VAX  ADA  COMPILER. 

Version  1.5  of  the  VAX  ada  compiler  was  used.  The  compiler  was  run  with 
full  optimization  with  regard  to  time.  This  optimization  includes  both 
local  and  global  optimizations  similar  to  those  which  are  performed  by 
the  MIPS  compilers  discussed  in  section  5.1.2. 

5. 2. 2. 2  VAX  C  Compiler 

The  VAX  C  compiler  compiler  can  perform  global  and  local  optimization 
by,  for  example,  doing  global  flow  analysis,  assigning  automatic 
variables  to  register  temporaries,  and  removing  invariant  computations 
from  loop,  to  mention  a  few.  The  compiler  also  does  peephole 
optimization.';  on  the  generated  machine  code. 

Version  2.3  of  the  VAX  C  compiler  was  used  in  this  study  and  run  with 
full  optimization  with  regard  to  time  (as  opposed  to  space). 
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6.  BENCHMARK  RESULTS/DISCUSSION. 


This  section  presents  the  benchmark,  results  obtained  by  running  both  the  Ada 
and  C  coded  FASTPNN  benchmark  discussed  in  section  3  using  the  coding  rules 
discussed  in  section  4  on  the  de^felopment  systems  discussed  in  section  5, 

It  is  important  to  reiterate  from  sections  5. 1.2.1  and  5. 1.2. 2,  that  the 
MIPS  C  code  was  compiled  and  executed  using  level  4  compiler  optimizations 
while  the  MIPS  Ada  code  was  compiled  with  level  1  optimizations.  Recall, 
that  this  is  due  to  the  fact  that  MIPS  level  3  and  4  compiler  optimizations 
were  not  supported  in  the  Ada  compiler  used  (Version  3.0).  Further,  when 
the  Ada  code  was  executed  using  level  2  optimizations  a  fatal  segmentation 
error  resulted.  Thus,  only  level  1  optimizations  could  be  used  for  the  Ada. 

In  an  attempt  to  circumvent  the  discrepancies  between  the  MIPS  C  compiler 
optimizations  and  the  MIPS  Ada  optimizations,  in  section  G.2,  projections  of 
Ada  execution  time  are  made  based  upon  level  4  compiler  assumptions.  The 
motivation  for  such  a  discussion  is  that  eventually  level  4  optimizations 
will  be  incorporated  in  Ada.  The  level  4  Ada  execution  time  projections  are 
based  upon  using  the  scaling  factor  associated  with  the  C  execution 
improvement  when  going  from  level  1  to  level  4  optimizations.  Thus,  the 
level  4  Ada  execution  time  projections  are  obtained  by  multiplying  the  Ada 
level  1  execution  time  by  the  C  level  4  execution  time,  and  then  dividing 
this  result  by  the  C  level  1  execution  time. 


6.1  Ada  Without  Checks  Versus  Ada  With  Checks 

To  support  the  use  of  exceptions,  Ada  performs  run  time  checks  to  determine 
whether  an  exception  should  be  raised.  In  practice,  a  clever  compiler  can 
determine  that  many  of  the  checks  can  be  safely  omitted.  Nevertheless,  a 
compiler  may  continue  to  generate  checks  that  a  programmer  knows  are 
unnecessary,  and  these  checks  may  make  a  critical  difference  in  the 
execution  time  of  a  program.  The  different  run  time  checks  include  access 
checks,  discriminant  checks,  index  checks,  length  checks,  range  checks, 
division  checks,  and  overflow  checks. 

To  investigate  the  execution  time  overhead  in  performing  checks,  the  FASTPNN 
Ada  code  was  executed  for  both  the  case  where  Ada  checks  were  performed  and 
also  the  case  where  all  Ada  checks  were  suppressed  (no  checks). 

Figure  8  contrasts  the  execution  efficiency  of  the  Ada  code  run  for  the  case 
where  Ada  run-time  checking  was  performed  versus  the  case  where  run-time 
checking  was  suppressed.  Figure  8  indicates  that  Ada  run-time  checks 
imposes  a  significant  penalty  on  the  execution  efficiency  of  the  Ada  for 
both  the  VAX  and  the  MIPS.  On  the  VAX,  the  PNN  execution  time  goes  from 
approximately  25  seconds  with  checks  off  to  approximately  36  seconds  with 
checks  on.  This  corresponds  to  a  43.'?  increase  in  time  for  performing  checks 
on  the  VAX.  On  the  MIPS,  the  execution  time  goes  from  approximately  2.4 
seconds  with  checks  suppressed  to  4  seconds  with  checks  preformed.  This 
corresponds  to  a  65%  increase  in  time  for  performing  checks  on  the  MIPS. 
Thus,  the  relative  penalty  for  performing  Ada  run-time  checks  is  moderately 
higher  for  the  HIPS  than  for  the  VAX. 
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6.2  C  Versus  Ada 


Figure  9  contrasts  the  execution  efficiency  of  the  C  code  with  the  Ada 
code.  The  relative  efficiency  of  C  versus  Ada  depends  on  which  machine 
the  code  is  executed  on.  If  the  code  is  executed  on  the  VAX,  than  the 
Ada  executes  more  efficiently  than  the  C.  For  the  case  of  the  VAX,  the 
Ada  executes  approximately  39. 5Z  quicker  than  the  C  when  Ada  checks  are 
suppressed  and  13.3%  quicker  than  the  C  when  Ada  checks  are  performed. 

On  the  other  hand,  if  the  code  is  executed  on  the  MIPS  the  C  code  runs 
quicker  than  the  Ada  code.  For  the  MIPS,  the  C  code  runs  approximately 
40%  quicker  than  the  Ada  code  when  Ada  checks  are  suppressed  and  114% 
quicker  than  the  Ada  code  when  Ada  checks  are  incorporated. 

Figure  10  contrasts  the  execution  efficiency  of  Ada  versus  C  for  the 
case  where  level  4  Ada  compiler  optimizations  are  projected  using  the 
scaling  technique  discussed  in  section  6.  Note  from  Figure  10  that  the 
Ada  execution  time  for  the  case  where  Ada  checks  are  not  incorporated  is 
nearly  identical  to  the  C  execution  tine. 

Based  on  Figures  9  and  10  it  is  concluded  that  there  is  little  or  no 
difference  between  the  inherit  execution  efficiency  of  Ada  (without  checks) 
and  that  of  C.  The  actual  execution  efficiency  of  Ada  versus  C  is  driven  by 
the  maturity  of  the  compilers  used  in  the  comparison.  With  regard  to  the 
MIPS,  the  fact  that  the  C  code  runs  quicker  than  the  Ada  is  attributed  to 
the  fact  that  the  present  MIPS  Ada  compiler  is  not  as  mature  as  the  MIPS  C 
compiler.  This  is  because  the  MIPS  Ada  compiler  was  run  with  level  1 
compiler  optimizations  while  the  MIPS  C  compiler  was  run  using  level  4 
compiler  optimizations.  Figure  10  indicates  that  eventually  when  the  MIPS 
Ada  compiler  matures  to  the  level  of  the  C  compiler,  the  Ada  code  will 
execute  in  nearly  identical  efficiency  (time)  as  the  C  code. 

The  fact  that  the  Ada  code  on  the  VAX  executes  more  efficiently  than  the 
C  ode  on  the  VAX  is  attributed  to  the  hypothesis  (this  was  not  proven 
:  this  study)  that  the  VAX  is  less  efficient  than  the  MIPS  in 

j.  forming  dynamic  allocation.  Thus,  in  the  case  of  Ada,  where  the 
»'  istrained  array  was  used  in  replacement  to  performing  dynamic  array 
aixocation  using  access  types,  the  resulting  Ada  code  was  more  efficient 
than  the  C  coae  where  dynamic  allocation  had  to  be  performed.  On  the 
MIPS,  'lovever,  where  dynamic  allocation  is  performed  more  efficiently 
this  .fference  between  Ada  and  C  was  nullified. 


6.3  MIPS  Versus  VAX 

I 

By  far  the  most  common  method  used  to  gauge  the  performance  of  a  particular 
machine  is  to  measure  its  performance  relative  to  the  VAX  11/780.  A  ratio, 
expressed  in  terms  of  VAX  MIPS,  is  obtained  by  dividing  the  time  required  to 
execute  a  given  benchmark  on  the  VAX  versus  the  time  to  execute  that  same 
benchmark  on  a  different  computer.  This  ratio  is  used  to  express  the 
machine's  performance  relative  to  the  VAX;  the  higher  the  ratio,  the  better 
the  machine's  performance.  Figure  11  contrasts  the  execution  efficiency  of 
the  MIPS  MAGNUM  3000  with  the  VAX  11/780.  Note  from  Figure  11  that  the 
ratio  of  MIPS  execution  time  versus  VAX  execution  time  is  highly  dependent 
on  whether  Ada  or  C  is  used.  If  C  is  the  language,  then  the  ratio  is 
approximately  22.2  VAX  MIPS.  If  Ada  is  the  language  the  MIPS  ratio 
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Figure  9.  FASTPNN  EXECUTION  TIME  COMIY.RISON 


Ada  (no  checks)  L-H  Ada  (checks) 


FIGURE  10.  FASTPNN  EXECUTION  TIME  COMPARISON  BASED 

UPON  ADA  LEVEL  4  OPTIMIZATION  PROJECTIONS 
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FIGURE  1 1 .  RELATIVE  THROUGHPUT  OF  MIPS  MAGNUM  300 

(NORMALIZED  TO  VAX  MIPS) 
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Ada  (no  checks)  Ada  (with  checks) 


drops  off  dramatically  and  is  dependent,  to  a  small  degree,  on  whether 
Ada  is  run  with  run-time  checks  or  without  run-time  checks.  If  Ada  is 
run  without  checks  than  the  ratio  is  approximately  10.3  VAX  MIPS.  In 
the  case  where  Ada  is  run  using  run-time  checks,  the  ratio  reduces  to 
approximately  9.0  VAX  MIPS. 

The  above  results  indicate  that  the  MIPS  MAGN'JM  3000  is  somewhere 
between  9  to  22.2  times  faster  than  the  VAX  11/780,  with  the  variance  in 
ratio  attributed  to  the  compiler.  Since  C  compilers  are  more  mature 
than  Ada  compilers,  it  is  concluded  that  the  C  run  time  results  are  a  better 
indicator  of  inherit  machine  performance  than  the  Ada  run  time  results. 

Thus,  it  is  concluded  that  the  actual  speedup  of  the  MIPS  over  the  VAX  is 
reflected  most  accurately  by  the  C  VAX  ratio  of  22.2  VAX  MIPS. 


6.4  MIPS  Profiler  Results 

The  purpose  in  profiling  is  to  help  identify  the  areas  of  code  where  most  of 
the  execution  time  -s  spent.  In  the  typical  program,  execution  time  is 
disproportionally  spent  in  relatively  few  sections  of  code.  Having 
identified  these  critical  sections  of  code,  it  is  profitable  to  improve 
coding  efficiency  in  those  sections. 

The  results  presented  in  this  section  were  obtained  by  using  the  MIPS  UNIX 
profiler.  In  both  the  C  and  Ada  case,  the  profiler  output  statistics  assume 
lOOX  cache  hits.  For  each  subroutine  in  PNN,  the  profiler  outputs  the 
following  statistics:  the  total  number  of  cycles  used  by  that  routine 
(CYCLES),  the  percentage  of  cycles  the  routine  uses  with  respect  to  the 
total  number  of  program  cycles  (^CYCLES),  the  cumulative  percentage  of 
cycles  (CUMX),  the  total  number  of  times  each  routine  is  called,  and  the 
number  of  cycles  used  by  the  routine  per  call  (CYCLES/CALL).  Note  that 
there  is  a  direct  relation  between  cycles  used  and  execution  time;  simply 
divide  the  cycles  used  by  the  clock  freqency  (25  MHz)  to  get  the  actual 
t  ime. 

Since  the  profiler  results  presented  in  this  section  are  categorized  by 
subroutine,  a  brief  description  of  the  seven  most  time  consuming  subroutines 
contained  in  FASTPNN  are  presented  in  Table  4.  Table  4  can  be  used  to  gain 
insight  into  the  time  breakdown  of  the  FASTPNN  routines  in  terms  of 
computational  functionallity . 

6.4.1  C  Profiler  Results 

In  this  section  the  profiler  was  run  on  the  FASTPNN  C  code  using  level  2 
compiler  optimizations.  Table  5  summarizes  the  profiler  results  for  the 
seven  most  time  consuming  subroutines  in  FASTPNN.  Table  5  lists  the  seven 
routines  in  descending  order  corresponding  to  their  overall  contribution 
with  regard  to  cycles  used  in  the  execution  of  FASTPNN.  Note  from  Table  5 
that  the  seven  routines  account  for  approximately  89X  of  the  total  cycles 
used  in  the  entire  FASTPNN  algorithm.  The  routine  which  accounts  for  the 
greatest  percentage  of  cycles  (time)  is  GetBucketStats .  From  Table  4  one 
can  see  that  GetBucketStats  is  concentrated  on  basic  floating  point 
mathematics  in  performing  mean  and  variance  calculations.  The  first  four  of 
the  routines  listed  in  Table  5  account  for  over  SOX  of  the  total  cycles 
used.  Also  note  from  Table  5  that  routine  Indxx  is  by  far  the  most  time- 
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TABLE  5.  MIPS  C  PROFILER  RESULTS 
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consuming  routine  with  regard  to  the  number  of  cycles  required  on  a  per  call 
basis,  requiring  over  970,000  cycles  per  call.  But  because  Indxx  is  only 
called  twice,  it  accounts  for  just  approximately  4.5%  of  the  total  cycles 
used. 

6.4.2  Ada  Profiler  Results 

In  this  section  the  profiler  was  run  on  the  Ada  code  using  level  1 
compiler  optimizations  and  suppressing  all  Ada  run  time  checks.  Table  6 
summarizes  the  profiler  results  for  the  seven  most  significant 
subroutines  in  FASTFNN.  Table  6  lists  the  seven  routines  in  descending 
order  corresponding  to  their  overall  contribution  with  regard  to  cycles 
used  in  the  execution  of  FASTPNN.  Note  from  Table  6  that  the  seven  routines 
account  for  approximately  94%  of  the  total  cycles  used  in  the  entire  PNN 
algorithm.  The  routine  which  accounts  for  the  greatest  percentage  of 
cycles  (time)  is  SplitBucket.  From  Table  4  one  can  see  that  SplitBucket 
is  computationally  intensive  in  the  areas  of  recursive  calling,  linked 
list  manipulation,  and  conditional  branching.  The  first  four  of  the 
routines  listed  in  Table  6  accounts  for  over  84%  of  the  total  cycles 
used  in  FASTPNN.  Also  note  from  Table  6  that  routine  Indxx  is  by  far 
the  most  time  consuming  routine  on  a  per  call  basis  requiring 
approximately  1.3  million  cycles  per  call. 

6.4.3  Ada  Versus  C  Profiler  Results 

Table  7  contrasts  the  performance  of  the  Ada  profiler  results  with  the  C 
profiler  results  for  each  of  the  main  routines  listed  in  the  previous 
two  sections.  The  seven  routines  are  presented  in  alphabetical  order  in 
Table  7.  Table  7  displays  the  amount  of  cycles  used  by  each  routine, 
the  relative  percentage  of  cycles  that  the  routine  contributes  to  the 
entire  FASTPNN  cycle  count,  and  a  ranking  of  each  routine  corresponding 
to  its  relative  contribution  of  cycles  used.  With  the  exception  of 
routine  Assesscandidate ,  the  results  displayed  in  Table  7  are  consistent 
from  the  standpoint  that  the  Ada  routines  use  more  cycles  than  the  C 
routines.  In  the  case  of  AssessCandidate,  the  C  requires  more  cycles 
than  the  Ada.  With  the  exception  of  routines  SplitBucket  and 
GetBucketStats,  the  relative  ranking  of  the  individual  routines  for  C 
versus  Ada  is  also  consistent.  In  the  case  of  routine  SplitBucket, 
SplitBucket  is  the  most  time  consuming  routine  (rank  1)  in  the  case  of 
Ada  while,  in  the  case  of  C,  SplitBucket  is  the  second  most  time 
consuming  routine  (rank  2).  In  the  case  of  routine  GetBucketStats, 
GetBucketStats  is  the  second  most  time  consuming  routine  (rank.  2)  in  the 
case  of  Ada,  while  GetBucketStats  is  the  most  time  consumirg  routine 
(rank  1)  in  the  case  of  C. 

The  last  row  of  Table  7  shows  the  total  cycles  from  the  seven  routines 
combined  and  the  cumulative  percentage  of  cycles  this  sum  comprises  of 
the  overall  cycles  used  in  the  FASTPNN  execution.  Note,  that  the  total 
cycles  used  by  the  Ada  exeeds  the  total  cycles  used  by  the  C  code.  This 
result  is  consistent  with  the  execution  results  previously  displayed  in 
Figure  9. 
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TABLE  6.  MIPS  ADA  PROFILER  RESULTS 
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TABLE  7.  COMPARISON  OF  MIPS  PROFILER  RESULTS 

(ADA  VERSUS  C) 
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7.0  PROJECTED  FASTPNN  REAL-TIME  REQUIREMENT 

The  goal  of  this  section  is  to  motivate  a  FASTPNN  real  time  requirement.  The 
basic  strategy  used  in  motivating  the  FASTPNN  real  time  requirement  derived  in 
this  section  is  outlined  below: 

A.  First,  determine  the  MBV  real  time  requirement/goal  (how  many 
targets  need  to  be  recognized  in  how  many  seconds?) 

B.  Next,  determine  what  percentage  of  time  an  MBV  system  will  be 
spend  on  FASTPNN 

C.  Last,  multiply  the  percentage  of  time  spent  on  FASTPNN  (B  above)  by 
the  overall  MBV  real-time  requirement  (A  above).  FASTPNN  must  then  be  able 
to  complete  its  execution  in  this  derived  interval  of  time 

There  is  currently  a  program  sponsored  by  the  Air  Force  called  Automatic 
Radar  Air-to-Ground  Acquisition  Program  (ARAGTAP)  which  is  focused  on 
establishing  a  real-time  MBV  capability.  The  ARAGTAP  goal  is  to 
identify  20  to  40  objects  (targets)  from  a  high  resoltion  SAR  image  in 
approximately  7  seconds.  If  we  conservatively  assume  a  20  to  30^  alarm 
rate,  this  requirement  translates  to  identifying  approximately  50  chips 
(of  which  20  to  40  may  be  actual  targets)  that  range  in  size  from  64  x 

64  pixels  to  128  x  128  pixels  in  7  seconds. 

Step  B  in  our  strategy  to  determine  a  FASTPNN  real-time  requirement  is  to 

determine  the  relative  percentage  of  the  time  that  the  MBV  algorithm 

will  spend  on  FASTPNN.  A  heuristic  reasoning  process  was  used  to  estimate 
that  FASTPNN  should  account  for  approximately  1.86X  of  the  total  MBV 
processing  time.  The  heuristic  reasoning  process  included  the  following 
assumptions: 

-  lOX  of  MBV  should  be  spent  on  prescreening/detection  and  the 
other  90%  should  be  spent  on  recognition 

-  of  the  remaining  90%  of  the  time  spent  on  recognition,  25%  of 
the  processing  time  should  be  spent  on  information  extraction  algorithms 
and  75%  should  be  spent  on  classification/matching  algorithms 

-  of  the  25%  of  the  time  spent  on  information  extraction, 

8.25%  of  this  time  should  be  spent  on  the  FASTPNN  algorithm. 

Mathematically  combining  all  of  the  above  assumptions  (by  multiplying), 
it  is  found  that  the  FASTPNN  algorithm  accounts  for  approximately  1.86%  of 
the  total  MBV  execution  time. 

Since  we  are  assuming  that  50  chips  must  be  processed  in  7  seconds,  it  is 
determined  that  the  FASTPNN  algorithm  must  be  able  to  process  a  single  chip 
(as  was  done  in  this  study)  in  .0026  seconds. 

In  section  6  the  best  case  FASTPNN  execution  time  was  1.87  seconds;  this  is 
nearly  3  orders  of  magnitude  slower  than  the  FASTPNN  real-time  requirement 
previously  derived  in  this  section. 


8.0  FASTPNN/PNN  IMPLEMENTATION  CONSIDERATIONS 


Based  upon  the  discussion  in  section  7,  it  is  evident  that  FASTPNN  execution 
efficiency  must  be  greatly  increased  if  it  is  to  function  within  real-time 
MBV  constraints. 

In  section  2.1  it  was  shown  that  FASTPNN  had  compuational  requirements 
that  are  O(NLOGN)  while  the  full  search  PNN  had  computational 
requirements  0(N**2).  The  fact  that  the  full  search  algorithm  is 
0(N**2)  and  that  FASTPNN  is  O(NLOGN)  does  not  necessarily  imply  that 
FASTPNN  will  execute  more  quickly  than  full  search  PNN  over  all  input 
data  sets.  In  fact,  for  up  to  relatively  large  data  set  sizes,  N,  it  is 
highly  conceivable  that  the  full  search  implementation  PNN  will  execute 
more  efficiently  than  FASTPNN.  This  is  attributed  to  the  fact  that 
there  are  initializations  and  other  overhead  associated  with  the  FASTPNN 
algorithm.  However,  as  N  gets  very  large,  the  volume  of  calculations 
associated  with  the  full  search  implementation  will  dominate  all 
overhead  associated  with  FASTPNN,  and  the  full  search  implemetation  will 
run  slower  than  FASTPNN.  A  logical  question  to  ask  is:  "for  a  given 
data  set  size  which  implementation  of  PNN  should  one  use"?  For  the  sake 
of  convenience  we  will  designate  the  parameter  BREAK_EVEN_N  to  refer  to 
the  data  set  size,  N,  where  the  execution  time  of  FASTPNN  and  full 
search  PNN  would  be  equal.  It  will  be  understood  that  for  N  < 

BREAK_EVEN_N ,  the  full  search  implementation  will  run  quicker  than 
FASTPNN,  and  for  N  greater  than  BREAK_EVEN_N ,  the  full  search 
implementation  will  run  slower  than  FASTPNN. 

It  is  important  to  emphasize  that  BREAK_EVEN_N  does  not  indicate  the 
data  set  size  where  both  implementations  of  PNN  are  equally  as  good. 

Given  equal  execution  time  for  both  implementations  of  PNN,  the  full 
search  implementation  of  PNIi  is  superior  to  FASTPNN  since  the  resulting 
quantization  vectors  are  optimum.  In  fact,  it  is  logical  to  assume  that 
significant  computational  savings  must  be  obtained  to  warrant  the  use  of 
FASTPNN  over  the  full  search  implementation. 

One  strong  disadvantage  of  the  full  search  PNN  is  the  large  memory 
requirement  which  would  be  required  if  a  large  input  data  set  size  is 
used.  For  example,  it  was  shown  by  example  in  section  2.1  that  if  the 
input  data  set  contain.^  1000  entries,  then  499,500  distance  calculations 
would  be  required  to  be  computed  in  the  first  iteration.  If  each 
distance  calculated  were  to  be  stored  in  one  word,  then  this  requirement 
translates  to  nearly  two  megabytes  of  memory. 

Intuitively,  the  parameter  BREAK_EVEN_N  is  dependent  on  the  actual 
architecture  used  in  executing  the  PNN  algorithm-  In  this  study,  where 
benchmarking  was  performed  on  general  purpose  data  processors  it  was 
assumed  (but  not  demonstrated)  that,  for  the  given  input  data  set, 

FASTPNN  was  more  time  efficent  than  full  search  PNN.  However,  if  one 
were  to  consider  executing  the  PNN  algorithm  on  a  signal  processor  (as 
opposed  to  a  data  processor),  the  choice  between  choosing  the  full  search 
implementation  of  PNN  versus  FASTPNN  could  certainly  be  altered.  A 
signal  processor  will  be  defined,  in  the  context  of  this  section,  as  a 
processor  specifically  designed  to  perform  sequences  of 
computations/operations  unaffected  by  actual  data  values.  Since  a  signal 
processor  is  more  suited  to  do  the  brute  force,  less  decision  intensive, 
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calculations  required  by  the  full  search  PNN  compared  to  the  more  data 
dependent  calculations  of  FASTPNN,  than  it  is  logical  to  conclude  that 
the  full  search  PNN  will  exhibit  higher  execution  efficiency  over  a 
larger  range  of  input  data  size. 

Figure  12  is  used  to  illustrate  the  points  made  in  the  above  paragraphs. 
Figure  12  shows  four  hypothetical  curves  which  are  used  to  contrast  the 
execution  efficiency  of  FASTPNN  versus  full  Search  PNN  for  both  a  data 

and  signal  processor  as  a  function  of  input  data  set  size,  N.  The  four 

curves  include  ;  (1)  the  execution  time  of  FASTPNN  on  a  signal  processor 
as  a  functjon  of  N;  (2)  the  execution  time  of  FASTPNN  on  a  data 
processor  as  a  function  of  N;  (3)  the  execution  of  FASTPNN  on  a  signal 

p’-ocessor  as  a  function  of  N,  and  (4)  the  execution  of  full  search  PNN 

cn  a  signal  processor  as  a  function  of  N.  Figure  10  shows  that  for 
small  values  of  N,  for  both  the  signal  and  data  processor 
implementation,  the  full  search  implementation  of  PNN  is  more  efficient 
than  FASTPNN.  But,  for  both  the  data  and  signal  processor  alike,  as  N 
gets  large,  FASTPNN  eventually  exhibits  higher  execution  efficiency  than 
the  full  search  implementation  of  PNN.  Note  from  Figure  12  that 
BREAK  EVEN  N  for  the  data  processor  occurs  for  a  much  smaller  input  data 
size,  than~it  does  for  the  signal  processor.  Also  note  from  Figure  12, 
that  the  data  processor  displays  slightly  better  execution  efficiency 
than  the  signal  processor  for  all  values  of  N  in  executing  FASTPNN.  But 
in  executing  the  full  search  implementation  of  PNN,  the  signal  processor 
significantly  outperforms  the  data  processor  for  all  values  of  N. 
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FIGURE  12.  HYPOTHETICAL  EFFICIENCY  COMPARISON  OF 

DATA  PROCESSOR  VERSUS  SIGNAL  PROCESSOR 
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9.  CONCLUSIONS 


9.1  Ada  With  Checks  Versus  Ada  Without  Checks 

Ada  run-time  checks  impose  a  significant  penalty  on  Ada  execution 
efficiency  for  both  the  MIPS  and  the  VAX.  On  the  VAX  there  was  a  43% 
relative  time  "penalty"  (increase  in  time)  associated  with  performing 
run-time  checks,  while  on  the  MIPS,  there  was  a  65%  time  penalty 
associated  with  performing  run-time  checks. 

9.2  C  Versus  Ada 

There  is  no  inherit  execution  efficiency  advantage  of  C  over  Ada  or  vice 
versa.  The  comparison  of  C  versus  Ada  depends  on  the  relative  maturity 
of  the  compilers  used.  For  instance,  on  the  MIPS  Magnum  3000,  the  C 
code  executes  more  efficiently  than  the  Ada  code.  While  on  the  VAX,  the 
Ada  code  executes  more  efficiently  than  the  C  code.  Thus,  the  resulting 
conclusion  is  that  for  the  HIPS  Magnum  3000,  the  C  compiler  is  more 
mature  than  the  Ada  compiler,  while  for  the  VAX  11/780  the  Ada  compiler 
is  more  mature  than  the  C  compiler. 


9.3  MIPS  Magnum  3000  Versus  VAX  11/780 


Depending  on  the  language  used  in  the  comparison,  the  MIPS  Magnum  3000 
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time)  than  the  VAX  11/780.  When  Ada  was  the  language  compared  on  both 
machines,  the  MIPS  executed  9  times  faster  than  the  VAX.  When  C  was  the 
language,  the  MIPS  executes  22  times  faster  that,  the  VAX.  The 
descrepancy  between  C  and  Ada  indicate  that  the  relative  efficiency  of 
MIPS  C  over  VAX  C  is  larger  than  the  efficiency  of  MIPS  Ada  over  VAX 
Ada. 


9.4  FASTPNN/  i  Real-Time  Implementation  Considerations 

The  best  case  FASTPNN  execution  time  obtained  from  executing  the  C  coded 
algorithm  on  the  MIPS  Magnum  3000  is  estimated  to  be  approximately  3 
orders  of  magnitude  too  slow  for  real  time  use.  Specialized  signal 
processor  hardware  will  be  required  to  boost  the  execution  efficiency 
of  FASTPNN  to  real  time  performance  levels.  When  using  specialized 
signal  processor  hardware,  it  may  be  advantageous  to  implement  the  full 
search  PNN  algorithm.  The  brute  force,  less  decision  intensive 
calculations  required  by  the  full  search  PNN  algorithm  make  it  more 
suitable  for  s^Tinal  essoi  hardware  application  than  the  FASTPNN 
algorithm.  ..ce  t'  ...oice  of  algorithm  implementation  is  dependent  on 
the  input  data  set  size,  a  study  should  be  performed  to  determine  which 
implementation  is  best  for  a  given  input  data  set  size.  Given  equal 
execution  times  for  both  the  full  search  PNN  implementation  and  the  FASTPNN 
implementation,  the  full  search  implementation  is  preferred  since  it 
provides  the  more  accu  results. 
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APPENDIX  A 


C  CODED  FASTPNN  BENCHMARK  SOURCE  CODE 


TABLE  OF  CONTENTS 


ROUTINE  PAGE 


MAIN  .  A-2 

FASTPNN  HEADER  FILE .  A-4 

FASTPNN  .  A-6 

VAX  TIMER  HEADER  FILE .  A-35 

VAX  TIMER  .  A-36 

MIPS  TIMER  HEADER  FILE .  A-37 

MIPS  TIMER .  A-38 

SAMPLE  OF  DATA  INPUT  FILE  .  A-39 


A-1 


ttinclude  "timer. h" 

#include  "fastpnn.h" 

^include  <stdio.h> 

float  *FastPNN  (float  *means,  float  ^weights,  iiit  count, 
int  dim,  int  ncntrds); 

Main(int  argc,  char  *-argv(  J) 

/* 

*  Function  name: 

*  main 

*  Purpose: 

*  Function  Main  is  the  driver  routine  for  the  Fast  Pairvise 

*  Nearest  Neighbor  Clustering  Algorithm  (FASTPNN).  This  routine 

*  calls  routines  to  perform  the  timing  of  the  FASTPNN  algorithm. 

*  This  routine  uses  command  line  arguments  to  pass  parameters  to 

*  itself  vhen  it  begins  executing.  At  the  command  line,  the  user 

*  enters  the  program  name  followed  by  the  following  parameters: 

*  the  input  file  name,  the  output  file  name,  and  the  number  of 

*  input  vectors. 

*/ 

{ 

float  *results; 
float  *positions; 
float  ^weights; 
float  a,b,c; 
int  count; 
int  dim; 

int  ncntrds,  i,  j,  k; 
int  loop_count  =  1; 

int  initial_dummy_time,  dummyarg,  dummy_elapsed_time; 

int  ini t_fastpnn_ time,  fastpnn_elapsed_time; 

float  fastpnn_iteration_time; 

int  vector_count ; 

int  weight  count ; 

int  position_count; 

FILE  *fpt_in; 

FILE  *fpt_out; 

i  =  0; 

j  =  0; 

fpt_in  =  fopGn(argv[l] ,  "r"); 
fpt_out  =  fopen(argv[2] ,  "w"); 
vectorcount  =  atoi(argv[3] ,  "r"); 
position_count  =  2  *  vector  count; 
ncntrds  =4;  “ 

positions  =  (float  *)calloc(position_count,  sizeof (float) ) ; 
weights  =  (float  *)calloc(vector_count ,  sizeof (float) ) ; 

while  (feof(fpt_in)  ==  0)  { 


A-2 


fscanf  (fpt_in,  "%f  %f  Xf",  &a,  &b,  &c); 
positions! 14^]  =  a; 
positionsi i+4-]  s  b; 
veightsi  j4-4-]  =  c; 


count  =  vGctor_count ; 

dim  =  positioncount/vector  count; 

printf( "enter  fas tpnnXn") ; 

ini tial_dummy_t ime  =  init_timer( ) ; 
printf ("got  initial_dummy_time\n") ; 
for  (i=0;  i<loop  count ;  i4-4-) 

{ 

dummy_arg  =  identi ty(dummy  arg); 

) 

dumrnyelapsedtime  =  elapsed_time(initial_dummy  time); 
printf ("dummy_elapsed_time  =  %d",  dummy_elapsGd_t ime) ; 

init_fastpnn_time  =  init_timer( ) ; 
for  (i=0;  i<loop  count;  i4-4-) 

{ 

results  =  FastPNN(positions,  weights,  count,  dim,  ncntrds); 
dummy_aig  =  identi  ty(duminy_arg); 

} 

fastpnn_elapsed_time  =  elapsed_time(init_fastpnn_time) ; 

fastpnn_iteration_time  =  (fastpnn_elapsed_time  -  dummy_elapsed_tiine)  / 

loop  count; 

fprintf ( fpt_out ,  "Fastpnn  executed  in  %f",  fastpnn_iteration  time); 
fprintf ( f pt_out ,  "microsecondsNn") ; 
f printf ( fptout ,  "Returned  resultsXn"); 
printarray(iesults,  dim,  ncntrds); 
f close( f ptout ) ; 


int  printarray( float  *array,  int  dim,  int  count) 


int  i,  j,  k; 

for  <i=0,  k=0;  i  !=  count;  i4-4-)  { 
printfC  ("); 

for  (j=0;  j  !=  dim;  j4-4-,  k4-4-)  { 
printfC'Xf",  array(k|); 
if  (j  !=  dim  -  1) 
printfC  "); 
else 

printf(")\n"); 

} 

) 

} 
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/*  This  is  the  header  file,  "Fastpnn.h",  which  contains  the 

*  data  structure  definitions,  function  prototype  definitions,  and 

*  symbolic  name  definitions  for  the  FastPNN  C  program  */ 


iifndef  _FastPNN 
ttdefine  _FastPNN 

#ifndef  TRUE 
^define  TRUE  -1 
#endif 

fifndef  FALSE 
♦define  FALSE  0 
♦endif 


♦define  KDNODE  0 
♦define  KDBUCKET  1 


♦define  KDMEMERR  "KD  tree  memory  allocation  error\n" 

♦define  BUCKETSIZE  8  /*  Number  of  entries  per  bucket  */ 
♦define  KDMERGE  0.5  /*  Fraction  of  buckets  merged  */ 
♦define  APTR  (char  *) 


struct  kdentry  { 


struct  kden 

int 

float 

float 

float 

float 

): 


ry  **next; 

splitlef t; 
weight ; 
*mean j 
*wmean ; 
*wsqmn; 


/*  k  dimensional  linked  list  pointers  */ 
/*  flag  used  for  bucket  splitting  */ 
/*  Weight  assigned  to  this  entry  */ 
/*  k  dimensional  sample  point  data  */ 
/*  k  dimensional  weighted  sample  data  */ 
/*  k  dim.  weighted  square  sample  data  */ 


struct  kdnode  { 


int 

dindx; 

/* 

Dimension  index 

*/ 

struct 

kdelem  *lower; 

/* 

Pointer 

to  kdeleras 

below 

thresh 

*/ 

struct 

kdelem  *upper; 

/* 

Pointer 

to  kdelems 

above 

thresh 

*/ 

); 


struct  kdbucket  { 
int  count; 
struct  kdentry  **lists; 
struct  kdentry  *entrya; 
struct  kdentry  *entryb; 
float  distort; 

}; 


/*  Cardinality  of  bucket  entries  */ 
/*  Pointers  to  sorted  data  linked  list  */ 
/*  First  element  of  candidate  pair  */ 
/*  Second  element  of  candidate  pair  */ 
/*  Distortion  induced  by  merging  pair  */ 


struct  kdelem  { 

int  type;  /*  value  of  KDNODE  or  KDBUCKET  */ 

union  {  /*  node  or  bucket  union  */ 

struct  kdnode  node; 
struct  kdbucket  bucket; 

]  norb; 

); 
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struct  kdtree  { 

*/ 
*/ 
■k/ 
*/ 


#endif 

void  fatal_message(char  string[]); 


lilt 

struct  kdelem 

int 

int 


dim;  /*  Dimension  of  tree  entries 

*root;  /*  Pointer  to  first  kd  tree  element 

nbuckets;  /*  Number  of  terminal  nodes 
nentries;  /*  Total  number  of  sample  points 
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^include  "fastpnn.h" 

float  *FastPNN  (means , weights , count , dim, ncntrds) 

float  *means,  ^weights; 
int  count,  dim,  ncntrds; 


/* 

*  Function  name: 

*  FastPNN 

* 

*  Purpose: 

*  Main  routine  for  the  Pairwise  Nearest  Neighbor  clustering  algorithm 

* 

*  Input  arguments; 


★ 

means 

-  sample  point  array 

■k 

weights 

-  sample  weight  array 

k 

count 

-  number  of  samples 

k 

dim 

-  dimensionality  of  the  sample  data 

k 

k 

ncntrds 

-  target  number  of  clusters  to  form  from  the  data 

*  Output  arguments: 

*  None 

* 


*  Returns: 

*  STACK  -  array  of  clusters 
*/ 


{ 

struct  kdentry  *entry; 

struct  kdtree  *tree,  *BuildKDtree(); 

float  *centroids; 

void  MergeDownKDtree( ) ,  DestroyKDtree() ; 

/*  char  *calloc();  */ 
int  i,  j,  k; 

if  (! (centroids  =  (float  *)calloc((unsigned)ncntrds*dim,sizeof (float)))) 
fatal_raessage(KDMEMERR) ; 

tree  »  BuildKDtree(means, weights, count, dim); 

MergeDownKDtree( tree, ncntrds) ; 

entry  =  tree->root->norb. bucket, lists(0] ; 
for  (i=0,k=0;  i!=ncntrds;  i++,entry=entry->next(01) 
for  (j=0;  j!=dim;  j++,k++) 

centroids[k]  =  entry->mean[ j J ; 

DestroyKDtree(tree) ; 

return(centroids) ; 

) 
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♦include  "fastpnn.h” 

struct  kdtree  *BuildKDtree(means, weights, count, dim) 

float  *means,  ^weights; 
int  count,  dim; 


/* 

*  Function  name: 

*  BuildKDtree 

■k 

*  Purpose: 

*  Constructs  an  initial  kD  tree  from  the  sample  data 

* 


*  Input  arguments: 


★ 

means 

-  sample  point  array 

★ 

weights 

-  sample  weight  array 

■k 

count 

-  number  of  samples 

★ 

dim 

-  dimensionality  of  the  sample  data 

*  Output  arguments: 

*  None 

* 

*  Returns: 

*  STACK  -  pointer  to  a  kD  tree 

*/ 


( 

struct  kdtree  *tree,  *CreateKDtree() ; 
struct  kdelem  *CreataFirstBucket() ; 

tree  =  CreateKDtree(dim) ; 

tree->root  =  CreateFirstBucket(means,weights,count,dim) ; 
tree->nbuckets  =  1; 
tree->nentries  =  count; 

return( tree) ; 

} 
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#include  "fastpnn.h" 
struct  kdtree  *CreateKDtree(dini) 
int  dim; 

/* 

*  Function  name: 

*  CreateKDtree 

* 

*  Purpose: 

*  Allocates  storage  for  a  kD  tree  data  structure 

* 

*  Input  arguments: 

*  dim  -  dimensionality  of  the  tree  data 

* 

*  Output  arguments: 

*  None 

* 

*  Returns: 

*  STACK  -  pointer  to  a  kD  tree 

*/ 

{ 

struct  kdtree  *treej 

if  (!(trs8  =  (struct  kdtree  *)ealloc((unsigned)l,sizeof(struct  kdtree)))) 
£atal_message(KDMEMERR) ; 

tree->dim  *  dim; 

return(tree) ; 

) 
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#include  "fastpnn.h" 
void  DestroyKDtree( tree) 


struct  kdtree  *tree; 

/* 

*  Function  name: 

*  DestroyKDtree 

* 

*  Purpose: 

*  Destroy  a  kD  tree 

* 

*  Input  arguments: 

*  tree  -  pointer  to  the  kD  tree 

* 

*  Output  arguments: 

*  None 

•k 

*  Returns: 

*  Nothing 
*/ 


{ 

void  DestroyLastBucketO; 
if  (‘ree->root) 

DestroyLastBucket(iree->root); 
cfree((char  *)tree); 


IS 


tinclude  "fastpnn.h" 

struct  kdelera  *CreateFirstBucket  (means,  weights,  count,  dim) 

float  *mean.s,  *veights; 
int  count,  dim; 


/* 

*  Function  name; 

*  CreateFirstBucket 

* 

*  Purpose: 

*  Create  and  initialize  the  first  bucket  in  a  kD  tree 

* 

*  Input  arguments: 


★ 

means 

-  sample  point  array 

★ 

weights 

-  sample  weight  array 

★ 

count 

-  number  of  samples 

★ 

dim 

-  dimensionality  of  the  sample  data 

★ 

*  Output  arguments: 

*  None 

* 

*  Returns; 

*  STACK  -  Initialized  bucket 
*/ 


{ 

struct  kdelem  *bucket,  *CreateKDbucket() ; 
struct  kdentry  *kdptr,  *lptr,  *CreateKDentry() ; 
void  SortBucketO; 
int  i,  j,  k; 

bucket  =  CreateKDbucket(dim) ; 

bucket->norb. bucket. count  =  count; 

for  (i=0,k=0;  i!=count;  i+i-)  { 

kdptr  =  CreateKDentry(dira) ; 

kdptr->veight  =  weightsli); 
for  <j=0;  jl=di'.;  j++,k++)  ( 
kdptr->mean[ j ]  =  meansfk]; 

kdptr->wmean[ j ]  =  kdptr->mean[ j )  *  kdptr->veight ; 
kdptr->vsqmn[ j ]  =  kdptr->mean[ j ]  *  kdptr->wmean[ j ] ; 

} 

if  (li)  { 

bucket->norb. bucket. lists[0]  =  kdptr; 

Iptr  =  kdptr; 

} 

else  { 

lptr->next [0]  =  kdptr; 

Iptr  =  kdptr; 

} 

} 

SortBucket  (bucket, dim) ; 
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reti  .n(buck,et) 


#lnclude  "fastpnn.h" 
struct  ltdelem  *CreateKDbucket(dim) 
int  dim; 


/* 

*  Function  name: 

*  CreateKDBuclcet 

* 

*  Purpose: 

*  Create  a  bucket  for  a  kD  tree 

■k 

*  Input  arguments: 

*  dim  -  dimensionality  of  the  sample  data 

k 

*  Output  arguments: 

*  None 

* 

*  Returns: 

*  STACK  -  Newly  created  bucket 

*/ 


{ 

struct  kdelem  ^bucket; 

/*  char  *calloc();  */ 

if  <! (bucket  =  (struct  kdelem  *)calloc( (unsigned)!, sizeof( struct  kdelem)))) 
fatal_message(KDMEHERR) ; 

if  ( I (bucket->norb. bucket .lists  »  (struct  kdentry  **) 
calloc( (unsigned)dim,sizeof(struct  kdentry  *)))) 
fatal_message(KDMEMERR) ; 

bucket->type  =  KDBUCKET; 

return(bucket) ; 

} 
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^include  "fastpnn.h" 
void  DesttoyLastBucket  (bucket) 
struct  kdelem  *bucket; 


/* 

*  Function  name: 

*  DestroyLastBuc.ket 

•k 

*  Purpose: 

*  Destroy  the  last  bucket  in  a  tree 

k 

*  Input  arguments: 

*  bucket  -  bucket  to  be  destroyed 

* 

*  Output  arguments: 

*  None 

k 

*  Returns: 

*  Nothing 
*/ 

{ 

void  DestroyKDbucket( ) ,  DestroyKDentry( ) ; 
struct  kdentry  *entry,  *next: 

entry  =  bucket->norb. bucket. lists(0J; 
while  (entry)  { 

next  =  entry->next(OJ ; 
DestroyKDen*;ry(entry) ; 
entry  =  next; 

} 

Des troyKDbucket( bucket ) ; 

) 
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tinclude  "fastpnn.h" 
void  DestroyKDbucket  (bucket) 
struct  kdelem  ^bucket; 


/* 

*  Function  name; 

*  DestroyKDBucket 

-k 

*  Purpose: 

*  Destroy  a  kD  tree  bucket 

* 

*  Input  arguments: 

*  bucket  -  bucket  to  be  destroyed 

* 

*  Output  arguments: 

*  None 

■k 

*  Returns; 

*  Nothing 
*/ 


{ 

cfree((char  *)bucket->norb. bucket .lists) 
cfree((char  *)bucket); 

} 
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#include  "fastpnn.h" 
struct  kdentry  •^CreateKDentry  (dim) 
int  dim; 


/* 

*  Function  name: 

*  CreateKDentry 

* 

*  Purpose: 

*  Create  a  kD  tree  bucket  entry  for  holding  a  sample  point 

* 

*  Input  arguments: 

*  dim  -  dimensionality  of  the  sample  data 

* 

*  Output  arguments: 

*  None 

* 

*  Returns: 

*  STACK  -  pointer  to  a  bucket  entry 

*/ 


struct  kdentry  *entry; 

/*  char  *calloc()j  */ 

if  (! (entry  >•  (struct  kdentry  *)calloc((unsigned)l,si2eof (struct  kdentry)))) 
fatal_message(KDM£MERR) ; 

if  ( ! (entry->next  =  (struct  kdentry  **) 

calloc((unsigned)diin,sizeof (struct  kdentry  *)))) 
fatal_message(KDMEMERR) ; 

if  ( ! (entry->mean  =  (float  *)calloc((unsigned)dim,sizeof(float)))) 
fatal_message(KDMEMERR) ; 

if  ( ! (entry->wmean  =  (float  *)calloc((unsigned)dim,sizeof(float)))) 
fatal_message(KDMEMERR) ; 

if  ( ! (entry->vsqmn  «  (float  *)calloc((unsigned)dim,si2eof(float)))) 
£atal_mes.sage(KDNEHERR) ; 

return(€ntry) ; 

) 
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•  include  "fastpnii.h" 
void  DestroyKDentry  (entry) 
struct  kdentry  *entry; 


/* 

*  Function  name: 

*  DestroyKDentry 

* 

*  Purpose: 

*  Destroy  a  kD  tree  bucket  entry 

* 

*  Input  arguments: 

*  entry  -  pointer  to  the  entry  to  be  destroyed 

* 

*  Output  arguments: 

*  None 

* 

*  Returns: 

*  Nothing 
*/ 


{ 

if  (entry)  { 

c£ree((char  *)entry->noxt) : 
cfree((char  *)entry->mean) ; 
c£ree((char  *)entry->wmean) ; 
cfree((char  *)entry->wsqmn) ; 
c£ree((char  *)entry); 

) 

) 
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#include  "fastpnn.h 


struct  kdeletn  *CreateKDnode( ) 

/* 

*  Function  name: 

*  CreateKDnode 

* 

*  Purpose: 

*  Create  a  k.D  tree  node 

* 

*  Input  arguments; 

*  Nothing 

* 

*  Output  arguments: 

*  None 

* 

*  Returns; 

*  STACK  -  pointer  to  the  newly  created  kD  tree  node 
*/ 


{ 

struct  kdelem  *node; 

/*  char  *calloc();  */ 

if  (!(riode  »  (struct  kdelem  *)calloc((unsigned)l,sizeof (struct  kdelem)))) 
fatal_message(KDMEMERR) ; 

node->type  =  KDNODE; 

return(node) ; 

) 
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♦include  "fastpnn.h" 
void  DestroyKDnode(node) 
struct  kdelem  *node; 


/* 

*  Function  name: 

*  DestroyKDnode 

* 

*  Purpose: 

*  Destroy  a  kD  tree  node 

*  Input  arguments: 

*  node  -  pointer  to  the  entry  to  be  destroyed 

* 

*  Output  arguments: 

*  None 

* 

*  Returns: 

*  Nothing 
*/ 


{ 

it  (node) 

cfree((char  ^^nodei: 

} 
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#inciude  "fastpnn.h" 

void  SortBucket  (elem,dim) 

struct  kdelera  *elem; 
int  dim; 


/* 

*  Function  name; 

*  SortBucket 

* 

*  Purpose: 

*  Sort  the  entries  in  a  kD  tree  bucket  across  each  dimension  separately 

* 

*  Input  arguments: 

*  elein  -  pointer  to  the  bucket  containing  the  entries  to  be  sorted 

*  dim  -  dimensionality  of  the  sample  data 

* 

*  Output  arguments: 

*  None 

-k 

*  Returns: 

*  Nothing 
*/ 


{ 

struct  kdentry  *entry,  **epbffr; 
void  indxxO; 

/*  char  *calioc();  */ 
int  i,  j,  count; 


count  =  elem->norb. bucket. count; 


if  (!(epbffr  =  (struct  kdentry  **)calloc((unsigne(?,)count , 

sizeof (struct  kdentry  *)))) 


fatal  message(KDMEMERR) ; 


for  (i=0;  i!=dim;  i++)  { 

for  (j»0,entry=elem->norb.bucket.lists[0] ;  jl-count;  j++)  { 
epbffr[ j ]  =  entry; 
entry  =  entry->next[OI ; 

] 

indxx(epbf fr ,&i) ; 

elem->norb. bucket. lists[il  »  epbffr[01; 
for  (j=l;  j<count;  j++) 

epbffr[j-l]->next(i]  =  epbffr[jj; 
epbffrf [count  -  lj]->next[i]  =  (struct  kdentry  *)0; 

} 

cfree((char  ^  pbffr); 

} 
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#include  "fastpnn.h" 

void  indxx(kdentry  **epbfrr,int  *i) 


/* 

*  Function  name: 

*  indxx 

if 

*  Purpose: 

*  Sort  an  array  of  indices  indx  based  on  the  data  arrin  using  an  indexed 

*  version  of  a  heap  sort.  Modified  from  Numerical  Recipes  indexx. 

if 

*  Input  arguments: 

*  arrin  -  array  of  data  used  for  sorting 

*  n  -  number  of  entries  in  arrin 

* 


*  Output  arguments: 

*  indx  -  Indices  specifiying  the  order  of  data  in  arrin 

* 

*  Returns: 

*  Nothing 
*/ 


{ 

int  1 J , ir, indxt, i ; 
float  q; 

for  (j=0; j<n;j++) 
indx[j]  =  jj 

1  =  n  »  1; 

ir  =  n  -  1; 

while  (TRUE)  { 
if  (1  >  0) 

C(  =  arrin[  (indxt=indx( — 1])]; 
else  { 

q  =  arrin( ( indxt=indx[ ir ] ) 1 ; 
indx[ir]=indxlO] 5 
if  (— ir  ==  0)  { 
indx{0]=indxt; 
return; 

} 

) 

i  =  1; 

j  ((1  +  1)^<  1)  -  1; 
while  (j  <=  ir)  { 

if  (j  <  ir  &&  arrin[ indx[ j 1 ]  <  arrin[ indx[ j+1 J ] ) 

j^+; 

if  (q  <  arrin[indx[j ]1)  { 
indx[ i]=indx[ j ] ; 
j  +=  ((i=j)  +  1); 

} 
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else  j=ir+l; 

} 

indx[i  l=iticixt ; 

} 

) 
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♦include  "fastpnn.h" 

void  MergeDownKDtree(tree,ncntrds) 

struct  kdtree  *tree; 
int  ncntrds; 


/* 

*  Function  name: 

*  MergeDownKDtree 

A" 

*  Purpose: 

*  Reduce  a  kD  tree  to  a  single  bucket  having  ncntrds  entries  using  the  PNN 

*  algorithm 

* 

*  Input  arguments: 

*  tree  -  pointer  to  a  kD  tree 

*  ncntrds  -  desired  number  of  entries  after  merging 

* 

*  Output  arguments: 

*  None 
■k 

*  Returns: 

*  Nothing 
*/ 


{ 

void  CompressKDtreeO  ,  BalanceKDtreeO ; 
struct  kdelem  *CollapseKI}node< ) ; 
int  nmerge,  ntile,  maxbkts; 
struct  kdelem  **bcktarray; 

/*  char  *calloc();  */ 

maxbkts  =  tree->nentries  /  (BUCKETSIZE  /  2); 


if  (!(bcktarray  =  (struct  kdelem  **)calloc( (unsigned)maxbkts , 

sizeof (struct  kdelem  *)))) 


£atal_message(KDMEMERR) ; 


while  (tree->nentries  >  ncntrds)  { 

BalanceKDt ree( tree , bcktarray ) ; 

ntile  =  (KDMERGE  *  tree->nbuckets  >  1)  ?  KDMERGE  *  tree->nbuckets  :  1; 
nmerge  =  ( ( tree->nentries  -  ncntrds)  <  ntile)  ? 

(tree->nentries  -  ncntrds)  :  ntllej 
CompressKDtree( tree, nmerge, bcktarray) ; 

} 

if  ( tree->root->type  ==  KDNODE) 

tree->root  =  Coi]apseKDnode( tree->root, tree->dim, &tree->nbuckets) ; 
cfree((char  *)bcktarray) ; 

} 
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tinclude  "fastpnn.h" 

void  BalanceKDtree( tree, bcktar ray ) 


struct  kdtree  *tree; 
struct  kdelem  **bcktarray; 


/* 

*  Function  name: 

*  RalanceKDtree 

* 

*  Purpose: 

*  Redistribute  the  entries  in  a  kD  tree  so  that  each  bucket  has 

*  approximately  the  same  number  of  entries 

•k 

*  Input  arguments: 

*  tree  -  pointer  to  a  kD  tree 

* 

*  Output  arguments: 

*  bcktarray  ~  array  used  to  retain  a  pointer  to  each  bucket  after  balancing 

k 

*  Returns: 

*  Nothing 
*/ 


( 

struct  kdelem  *CollapsekDnode( ) ,  *SplitBucket() ,  **bptr; 
float  *mean,  *vvar; 

/*  char  *calloc();  */ 

if  (!(mean=  (float  *)calloc((unsigned)tree->dim,sizeof (float)))) 
fatal_message(KDHEMERR) ; 

if  (!(vvar  =  (float  *)cailoc( (unsigned) tree->dim, sizeof( float ))) ) 
fatal_message(KDMEMERR) ; 

if  ( tree->root->type  ==  KDNODE) 

tree->root  =  CollapseKDnode( tree->root, tree->dim,&tree->nbuckets) ; 

bcktarray [0]  =  tree->root) 
bptr  M  &bcktarray(li ; 

tree->root  = 

SplitBucket( tree->root, tree->dim,&tree->nbuckets,iibptr ,mean, wvar) ; 

cfree((char  *)mean); 
cfree((char  *)vvar); 
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tinclude  "fastpnn.h" 

struct  kdelera  *CollapseKDnode(elem,dim, bctr) 

struct  kdelem  *elem; 
int  dim,  *bctr; 


/* 

*  Function  name: 

*  CollapseKDnode 

* 

*  Purpose: 

*  Recursive  function  used  to  collapse  the  entries  in  a  kD  tree  into  a 

*  single  bucket 

•k 

*  Input  arguments: 

*  elem  -  pointer  to  a  kD  tree  node  to  be  collapsed 

*  dim  -  dimension  of  the  data  within  the  tree 

* 

*  Output  arguments: 

*  bctr  -  pointer  to  a  counter  used  to  keep  track  of  the  total  number 

*  of  buckets  in  the  tree 

* 

*  Returns: 

*  STACK  -  pointer  to  the  bucket  resulting  from  the  collapse 

*/ 


( 

struct  kdentry  **entptr,  *lentry,  *rentry; 
struct  kdelem  *bucket; 
void  DestroyKDnode( ) ,  DestroyKDbucket( ) ; 
int  i; 

if  (elem->norb. node. lower-> type  !=  KDBUCKET) 

elem->norb. node. lover  =  CollapseKDnode(elem->norb. node. lover, dim, bctr) ; 
if  (elem->norb. node. upper-> type  !=  KDBUCKET) 

elem->norb. node. upper  =  CollapseKDnode(elem->norb.node.upper,dim,bctr) ; 

bucket  =  elem->norb. node. lover; 

for  (i»0;  il=dlm;  i++)  { 

entptr  =  &(elem->norb. node. lower->norb. bucket . llsts( i) ) ; 
lentry  =  elem->norb. node. lower->norb. bucket . lists[ i j ; 
rentry  =  elem->norb.node.«pper->norb. bucket . lists[ i] ; 
while  (lentry  &&  rentry)  ( 

if  (lentry->mean[ i ]  <  rentry->mean[ i ] )  { 

*entptr  =  lentry; 
entptr  =  6i(lentry->next[i]); 
lentry  =  lentry->next[ii ; 

) 

else  ( 

*entptr  =  rentry; 
entptr  =  &(rentry->next[i]); 
rentry  =  rentry->next [ i ] ; 
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} 

if  (lentry)  { 
while  (lentry) 


{ 


*entptr  =  lentry; 
entptr  =  &(lentry->next ( i j ) ; 
lentry  =  lentry->next[i] ; 

} 

} 

if  (rentry)  { 
while  (rentry)  { 

*etitptr  =  rentry; 
entptr  =  &( rentry->next ( i ] ) ; 
rentry  =  rentry->next [ i ) ; 

} 

1 


) 

bucket->norb. bucket . count  =  elem->norb. node. lower->norb. bucket . count  + 

elem->norb.no<le.upper->norb.  bucket  .count ; 
Dest royKDbucket( elem->norb. node .upper ) ; 

DestroyKDnode  (elem) ; 

(*bctr)  — ; 
return(bucket) ; 
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#include  "fastpnn.h" 

struct  kdelem  *SplitBuck.et(oldbucket, dim, bctr,bptr, mean, war) 


struct  kdelem  *oldbucket,  ***bptr; 
int  dim,  *bctr; 
float  *mean,  *wvar; 


/* 

*  Function  name: 

*  SplitBucket 

■k 

*  Purpose; 

*  Recursive  function  used  to  split  a  kD  tree  bucket  into  two  smaller 

*  buckets  having  half  as  many  entries 

* 

*  Input  arguments: 

*  oldbucket  -  pointer  to  kD  tree  bucket  to  be  split 

*  dim  -  dimension  of  the  data  within  the  tree 

*  mean  -  scratch  array  used  for  calculating  bucket  means 

*  wvar  -  scratch  array  used  for  calculating  bucket  weighted  variances 

* 


*  Output  arguments: 

*  bctr  -  pointer  to  a  counter  used  to  keep  track  of  the  total  number 

*  of  buckets  in  the  tree 

*  bptr  -  pointer  to  an  array  of  pointers  to  buckets  in  the  tree 

* 


*  Returns: 

*  STACK  -  pointer  to  the  node  resulting  from  the  split,  or  the  original 

*  bucKet  if  the  number  of  entries  in  the  bucket  is  small  enough 
*/ 


{ 

struct  kdelem  *newnode,  *newbuckct,  *CreateKDbucket( ) ; 
struct  kdentry  **oldptr,  **newptr,  *entry; 
int  i,  j,  bcount,  medindx; 
void  GetBucketStats( ) ; 

if  (oldbucket->norb. bucket. count  >  BUCKET5IZE)  { 

Ge tBucketStats( oldbucket , dim, mean, war); 

for  (i=l,j=0;  i<dim;  i++) 
if  (vvar(i]  >  war[j]) 

j  «  i; 

bcount  =  oldbucket->norb. bucket. count; 

medindx  =  (bcount  +  1)  /  2;  /*  Uneven  splits  go  left  */ 

newnode  =  CreateKDnode() ; 

newnode->norb.node.dindx  ^  j; 

newbucket  =  CreateKDbucket(dim) ; 

newnode->norb. node . lower  =  oldbucket; 
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nevnode->norb. node. upper  *  newbucket; 

(  ’^bctr)  ++; 

(**bptr)  =  newbucket; 

(  *bptr)  ++; 

for  (1=0,  entry=oldbucket->norb. bucket. listsi j ) ;  i<medindx;  i++)  { 
entry->splitlef t  =  TRUE; 

entry  =  entty->next ( j ) ; 

} 

for  (i=medindx;  i<bcount;  i++)  { 
entry->splitlef t  =  FALSE; 

entry  =  entry->next[ j 1 ; 

} 

oldbucket->norb. bucket . count  =  medindx; 
newbucket->norb. bucket . count  =  bcount  -  medindx; 

for  (1=0;  i!=dim;  i++)  { 

oldptr  =  &oldbucket->norb. bucket. lists[i] ; 
newptr  ^  iinevbucket-->norb. bucket. lists[i] ; 
entry  =  oldbucket->norb. bucket. listsli] ; 
while  (entry)  { 

if  (entry->splitlef t)  { 

*oldptr  =  entry; 

oldptr  =  (struct  kdentrv  **)(&entty->next I i j ) : 

) 

else  { 

*newptr  =  entry; 

newptr  =  (struct  kdentry  **)(&entry->next[i] ) ; 

) 

entry  =  entry->next[i] ; 

} 

*oldptr  =  (struct  kdentry  *)0; 

*newptr  =  (struct  kdentry  *)0; 

} 

newnode->norb. node. lower  » 

SplitBucket(newnode->norb.node. lower, dim, bctr, bptr, mean, wvar) ; 
newnode->norb. node. upper  » 

Spli tBucket (newnode->norb. node. upper, dim, bctr , bptr, mean, wvar) ; 

return(nevnode) ; 

} 

else 

return(oldbucket) ; 
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tinclude  "fastpnn.h" 

void  GetBucketStats  (elem, dim, mean, war) 

struct  kdelem  *elem; 
int  dim; 

float  *mean,  ^wvar; 


/* 

*  Function  name: 

*  GetBucketStats 

* 

*  Purpose: 

*  Calculate  the  k  dimensional  means  and  weighted  variances  for  a  bucket 

* 

*  Input  arguments: 

*  elem  -  pointer  to  the  kD  tree  bucket  for  which  the  statistics  are 

*  to  be  calculated 

*  dim  -  dimension  of  the  data  within  the  tree 

* 


*  Output  arguments: 

*  mean  -  array  used  for  calculating  bucket  means 

*  wvar  ~  array  used  for  calculating  bucket  weighted  variances 

•k 

*  Returns: 

*  Nothing 

*  / 


{ 

struct  kdentry  *entry; 
float  wgtsum; 
int  i; 

wgtsum  =  0.0; 

for  (i=0;  i!=dim;  i++)  { 
mean[i]  =  0.0; 
wvar [ i )  =  0.0; 

} 

entry  =  elem->norb. bucket. lists [0 1 ; 
while  (entry)  { 

wgtsum  +=  entry->weight ; 
for  (i=0;  i!=dira;  i++)  { 
mean[i]  +=  entry->vmean[ i ] ; 
vvar[i]  +=  entry->vsqmn[i] ; 

} 

entry  =  entry->next[0] ; 


for  (i=0;  i!=dim;  i++)  { 
mean[i]  /=  wgtsum; 

wvar[i]  =  (vvar[i]  /  wgtsum)  -  (raeanli]  *  meanlij); 

} 
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#include  "fastpnn.h" 
static  int  PucketConipare(one,  two) 
struct  kdelem  **one,  **tvo; 


/* 

*  Function  name: 

*  BucketCompare 

* 


*  Purpose: 

*  Function  used  by  the  UNIX  qsort  routine  to  compare  merge  distortions  of 

*  two  buckets 

* 

*  Input  arguments: 

*  one  -  pointer  to  a  pointer  to  the  first  kD  tree  bucket  used  in 

*  the  comparison 

*  two  -  pointer  to  a  pointer  to  the  second  kD  tree  bucket  used  in 

*  the  comparison 


*  Output  arguments: 

*  None 

* 


*  Returns: 

*  STACK 

* 

*/ 


the  value  1 
the  value  0 
the  value  -1 


if  di.stortion(one)  >  distortion(tyo) 
if  distortion(one)  =  distortionf two) 
if  distortion(one)  distortion( two) 


if  ( (*one)->norb, bucket .distort  <  (*two)->norb. bucket .distort) 

return(-l) ; 

else  it  ( (*one)->norb. bucket .distort  ==  (*two)->norb. bucket .distort) 
return(  0); 
else 

return(  1); 
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#include  "fascpnn-h" 

void  CompressK-Dtree  (tree,  nmecge,  bcktarray) 

struct  kdtree  *tree; 
int  nmerge; 

struct  kdelem  **bcktarray; 


/* 

*  Function  name: 

*  CompressKDtree 

* 

*  Purpose: 

*  Function  used  to  merge  bucket  entry  pairs  into  single  bucket  entries  for 

*  a  fixed  fraction  of  the  total  number  of  buckets 

•k 

*  Input  arguments: 

*  tree  -  pointer  to  the  kD  tree  undergoing  the  merge 

*  nmerge  -  number  of  bucket  pairs  to  merge 

*  bcktarray  -  array  of  pointers  to  all  buckets  in  the  tree 

* 

*  Output  arguments: 

*  None 

* 

*  Returns: 

*  Nothing 

* 

*/ 


{ 

void  AssessCandidateO ,  ReduceKDbucket() ; 
int  i,  ncount,  BucketCompareC ) ; 

for  ( i=0,ncount=0;  i ! =tree->nbuckets;  i++) 

AssessCandidate( bcktarray!  i  ] ,  tree->diiii,  bcktarray ,  incount ) ; 

/*  Handle  end  game  situations  */ 

if  (nmerge  >  ncount) 
nmerge  =  ncount; 

if  (ncount  >  1) 

qsort((char  *)bcktarray, ncount ,sizeof (struct  kdelem  *) , BucketCompare) ; 

for  (i=0;  i!=nmerge;  i++) 

ReduceKDbucket(bcktarray{ i ] , tree->dim) ; 

tree->nentries  -=  nmerge; 
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#include  "fastpnn.h" 

void  Asses. sCandidate(el.eiTi,  dim,  bar  ray  ,ncount) 

struct  kdelem  *eleiu,  **barray; 
int  dim,  *ncount; 


/* 

*  Function  name: 

*  AssessCandidate 

* 

*  Purpose: 

*  Determine  the  minimal  distortion  than  can  be  produced  by  merging  a  pair 

*  of  bucket  entries 

* 

*  Input  arguments; 

*  elem  -  pointer  to  the  kD  tree  bucket  under  evaluation 

*  dim  -  dimension  of  the  data  within  the  entries 

* 

*  Output  arguments: 

*  barray  -  array  of  pointers  to  all  buckets  that  can  be  merged 

*  ncount  -  pointer  to  counter  used  to  keep  track  of  the  total  number 

*  of  buckets  that  can  be  merged 

* 

*  Returns: 

*  Nothing 

* 

*/ 


struct  kdentry  *ientry,  *jentry; 
float  reduction,  dotprd,  diff; 
int  i,  j,  k,  firsttime; 

first  time  =  TRUE; 

if  (elem->norb. bucket. count  >  1)  { 

ientry  =  elem- >norb. bucket . 1 ists[0] ; 

for  (i=0;  i<elem->norb. bucket. count-1;  i++, ientry=ientry->next [0] )  { 
jentry  =  ientry->next [0] ; 

for  j<elem->norb. bucket  .count;  j++,  jentry*:jentry->next[0])  { 

for  (k=0, dotprd=0.0;  k<dim;  k++)  { 

diff  =  ientry->mean[kj  -  jentry->mean(k] ; 
dotprd  +=  diff  *  diff; 

1 

reduction  =  dotprd  *  ientry->weight  *  jentry->weight  / 

(ientry ->veight  +  jentry->weight) , 
if  ((reduction  <  elem->norb. bucket. distort)  | |  firsttime)  { 
elem->norb. bucket. distort  «=  reduction; 
elem->norb. bucket . enttya  =  ientry; 
elem->norb. bucket . entryb  =  jentry; 
firsttime  =  FALSE; 

) 

) 

} 

barrayl (*ncount)++]  =  elem; 
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tinclude  "fastpnn.h” 

void  ReduceKDbuck.et(elem,dim) 

struct  kdelem  *elem; 
int  dim; 


/* 

*  Function  names 

*  ReduceKDbucket 

* 

*  Purpose: 

*  Merges  a  pair  of  bucket  entries  into  a  single  entry 

•k 

*  Input  arguments: 

*  elem  -  pointer  to  the  kD  tree  bucket  whose  entries  are  to  be  merged 

dim  -  dimension  of  the  data  within  the  entries 

* 

*  Output  arguments: 

*  None 

k 

*  Returns: 

*  Nothing 

* 


{ 

struct  kdentry  *ientry,  *jentry,  **leptr; 
float  newveight; 
int  i,  rmcnt; 

ientry  =  elem->norb. bucket. entrya; 
jentry  =  elera->norb. bucket. entryb; 

/*  Remove  ientry,  jentry  from  the  list  */ 

cor  (i=0;  iUdim;  i+-t-)  { 

for  (leptr  =  &(elem->norb. bucket . lists! i ]), rmcnt=0;  rmcntI-2;)  { 
if  ((*leptr==»jentry)  ||  (*leptr==ientry) )  { 

*leptr  =  (*leptr)->ne:ct[  i  J ; 
rmcnt++; 

} 

else 

leper  =  &( (*leptr )->next [ i] ) ; 

) 

} 

iieweight  =  ientry->weight  +  jentry->weight ; 
for  (i=0;  i!=dim;  i+-(-)  { 

ientry->mean[ i 1  =  ( ientry->mean[ i ]  *  ientry->veight  + 

jentry->mean(i)  *  jentry->weight )  /  newveight; 
ientry->vmean[i  ]  =  ientry->i,iean(  i  ]  *  newveight; 
ieiitry->wsqmn( i 1  =  ientry->mean| i ]  *  ientry->vmean[ i ] ; 

} 

ientry->weight  =  newveight; 
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/*  Reinsert  ientry  into  the  list  in  the  proper  order  */ 

for  (i=0;  i!=dim;  i++)  { 

for  (leptr  =  &(elem->norb. bucket. listsli]); 

(*leptr  &&  ((*leptr)->tnean(i]  <  ientry.->mean[  i  ] ) )  ; 
leptr  =  &((*leptr)->nextli| )); 
ientry->next[i]  =  *leptr; 

*leptr  =  ientry; 

} 

eleni->norb. bucket . count-- ;  /*  Decrement  the  total  entry  count  */ 

DestroyKDentry(j entry);  /*  Free  up  jentry  memory  */ 


void  fatal  message(char  string[]) 

{ 

print f("Xs\n" .string) ; 
return; 

} 
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/*  This  is  the  header  file, "timer .h" ,  which  contains  the  declarations 

*  (or  function  prototypes)  for  the  VAX/VMS  timing  routines 

* 

*/ 

int  init_timer(void) ; 

int  elapsed_time<int  starting_time) ; 

int  identity  (int  arg); 
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#include  "timer. h" 
ttinclude  <time.h> 

#include  <stdio.h> 

/*  This  file  contains  the  VAX/VMS  timing  routines.  Note  that  the 

*  functions  init_timer  and  elapsed_time  make  use  of  a  predefined 

*  C  function  called  "times".  Function  times  returns  the  accumulated 

*  CPU  time  in  a  predefined  time  structure  called  tbuffer_t. 

•k 

*/ 

tbuffer_t  *init_time_ptr,  *f inal_time_ptr; 
tbuffer_t  buffer,  init_time,  final_time; 

int  init_timer  (void) 

{ 

int  current_time; 

ini t_time_ptr  =  &init_time; 

times( ini t_time_ptr) ; 

current_time  =  (init_time_ptr  ->  proc_user_time)  *  10000; 
return  current_time; 

} 


int  elapsed_time(int  start_tirae) 

{ 

int  time_elapsed; 
tinai_time_ptr  =  &final_time; 
times(f  inal_tiiiie_ptr) ; 

time_elapsed  =  ((final  time_ptr->proc  user  time)  *  10000)  -  start_time 
return  time_elapsed;  ~  ” 

) 

int  identity  (int  arg) 

{ 

int  some_value  =  0; 
some_value  =  some  value  +  arg; 
return  somevalue; 

} 


A-36 


/*  This  is  the  header  file,  "timer. h",  vhich  contains  the  declarations 
*  (or  function  prototypes)  for  the  MIPS/UNIX  timing  routines. 

*/ 


long  int  init_ti!ner(void) ; 
long  int  elapsed_time(void) ; 
int  identity  (int  arg) ; 
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#lnclude  "timer.h'' 
tinclude  <time.h> 

#include  <stdio.h> 

/*  This  file  contains  the  MIPS/UNIX  timing  routines.  Note  that  the 

*  functions  inittimer  and  elapsed_time  make  use  of  the  predefined 

*  C  function  called  "clock”.  Function  clock  returns  the  amount  of 

*  CPU  time  used  since  the  first  call  to  clock. 

* 

*/ 

long  int  init_timer  (void) 

{ 

long  int  current_time; 
current  time  =  cTock() ; 
return  current  time; 

} 


long  int  elapsed_time(void) 

{ 

long  int  timeelapsed; 
time_elapsed  =  clock() ; 
return  time  elapsed; 


int  identity  (int  arg) 

I 

int  some_value  =  0; 
some_value  «  some_value  -t-  arg; 
return  some_value; 

} 
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i 


58.000000  19.000000  62.000000 
59.000000  19.000000  74.000000 
60.000000  19.000000  70.000000 
57.000000  20.000000  76.000000 
58.000000  20.000000  90.000000 
59.000000  20.000000  91.000000 
60.000000  20.000000  78.u00000 
61.000000  20.000000  52.000000 
56.000000  21.000000  75.000000 
57.000000  21.000000  96.000000 
58.000000  21.000000  101.000000 
59.000000  21.000000  91.000000 
60.000000  21.000000  67.000000 
56.000000  22.000000  75.000000 
57.000000  22.000000  100.000000 
58.000000  22.000000  108.000000 
59.000000  22.000000  104.000000 
60.000000  22,000000  92.000000 
61.000000  22.000000  72.000000 
56.000000  23.000000  84.000000 
57.000000  23.000000  135.000000 
58.000000  23.000000  159.000000 
59.000000  23.000000  165.000000 
60.000000  23.000000  155.000000 
61.000000  23.000000  124.000000 
62.000000  23.000000  61.000000 
56.000000  24.000000  105.000000 
57.000000  24.000000  164.000000 
58.000000  24.000000  190.000000 
59.000000  24.000000  196.000000 
60.000000  24.000000  185.000000 
61.000000  24.000000  151.000000 
62.000000  24.000000  74.000000 
55.000000  25.000000  59.000000 
56.000000  25.000000  109.000000 
57.000000  25.000000  173.000000 
58.000000  25.000000  200.000000 
59.000000  25.000000  206.000000 
60.000000  25.000000  194.000000 
61.000000  25.000000  160.000000 
62.000000  25.000000  81.000000 
63.000000  25.000000  58.000000 
36.000000  26.000000  52.000000 
56.000000  26.000000  105.000000 
57.000000  26,000000  165.000000 
58.000000  26,000000  191.000000 
59.000000  26.000000  197.000000 
60.000000  26.000000  185.000000 
61.000000  26.000000  152.000000 
62.000000  26,000000  80.000000 
63.000000  26.000000  58.000000 
36.000000  27.000000  60.000000 
37.000000  27.000000  73.000000 
40.000000  27.000000  85.000000 
41.000000  27.000000  94.000000 
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vith  timer; 

vith  FASTPNN; 

vith  READDATa; 

vith  text_io;  use  text_io; 

vith  Data  Struct_Pkg;  use  Data_Struct_Pkg; 

procedure  main  is 


—  Procedure  Name:  Main 

—  Purpose  :  This  is  the  main  routine  for  the  FASTPNN  algorithm. 

The  user  is  prompted  to  input  the  name  o£  the  input 
file,  the  name  of  the  output  file,  and  the  number  of 
input  vector  data  values.  External  procedure 
READDATA  is  called  by  Main  to  read  in  the  values 
from  the  input  data  file.  Main  performs  the  timing 
of  the  FasiPNN  algorithm  (by  calling  timing 
routines  contained  vithin  Timer_Picg)  and  outputs  the 
results. 


LFUGTHIN  :  INTEGER; 

LENGTHOUT  :  INTEGER; 

POSITIONSLAST  :  INTEGER; 

RESULTS_LAST  i  INTEGER; 

WEIGHTS_LAST  :  INTEGER; 
count  !”integer; 
dim  :  integer; 
ncntrds  :  Integer; 

INNAME  :  STRING(1  ..  80); 

OUTNAHE  :  STRING (1  ..  80); 

VECT0R_NUM  ;  INTEGER; 

loopccunt  :  constant  1; 
fastpnn_elapsed  :  integer; 

£astpnn_timer  :  timer .microsec_timer; 

dummy  timer  ;  timer .microsec_timer; 
dummy _elapsed_time  :  integer; 
dummy"arg  :  integer; 

INFILE  :  FILETYPE; 

OUTFILE  s  FILE_TYPE; 

package  float_io  is  nev  text_lo.float_io( float); 
package  int_lo  is  nev  text_io.integer_io( integer); 

procedure  printarray  (OUTFILE  :  in  out  FILE_TYPE; 

arrayy  ;  in  means  array_type; 

dim  :  in  integer;~count  :  in  integer)  is 


k  :  integer; 

package  float_io  is  nev  text_io.float_io( float); 

begin 

k  :*  0; 
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for  i  in  0  . .  count  -  1  loop 
put<OUTFILE,  "  ("); 

for  j  in  0  . .  dim  -  1  loop 

float_io.put(OUTFILE,  arrayy(k,)); 
k  k,  +  1; 
if  (j  /-  dim-1)  then 
put<OUTFILE,  "  •'); 

else 

put_line(OUTFILE, 
end  if; 
end  loop; 
end  loop; 
end  printarray; 

begin  —  MAIN 

PUT_LINE( "Enter  Name  of  Input  File.  "); 

GET_LINE( INNAME ,  LENGTHIN) ; 

PUT  LINE( "Enter  Name  of  Output  File  "); 

GET"LINE(0UTNAME,  LENGTHOUT); 

PUT_LINE< "Enter  the  number  of  input  vectors  "); 

INT_IO . GET ( VECTOR_NUM ) ; 

OPENdNFILE,  IN^FILE,  INNAME(1  ..  LENGTHIN)); 

CREATE(OUTFILE,  0UT_FILE,  0UTNAME(1  ..  LENGTHOUT)); 

POSITIONS  LAST  2  *  VECTOR  NUM  -  1; 

RESULTS  LAST  :*  2  *  VECTOR  NUM  -  1; 

VEIGHTS'LAST  V£CT0R_NUM''-  1; 

DECLARE 

positions  :  means_array_typc(0  ..  POSITIONS_LAST) ; 
results  5  raeans_array_type(0  ..  RESULTS_LAST) ; 
weights  :  veights_array_type(0  WEIGHTS_LAST) ; 
begin 

count  ;»  weights' length; 

dim  im  positions'length/count; 

READ_DATA  (INFILE,  positions,  weights); 

ncntrds  :=«  4;  —  number  of  desired  output  vectors 

—  first  time  the  dummy  loop 

timer. ini t_timer(dummy_timer) ; 
for  i  ir  1  . .  loop_count  loop 

dummy  arg  timer . identi ty(duinmy_arg) ; 
end  loop; 

if  timer. always_true  then 

dummy _elapsed_tlme  :■  timer .elapsed  tlme(dummy  timer) 
end  if;  ~  ~ 
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—  now  time  the  toulinc  ol  interest 


timer. init_tiBei( fas tpnnt leer ), 
for  i  in  1  . .  loopcount  loop 

Fas  t  [’NN(  pos  i  i  ions  ,  weights,  cownt.  di*.  ncnttds.  resuiti); 
duauay  arg  :«  t  i»er  ■  ioent  i  ty((iuM>  arg  ' ; 
end  loop; 

if  t iaet . always  t rue  then 

f  as  tpnnclapseo  :•  t  t«e:  .elapsrc  t  las  :paA_  1 > . 
end  It; 


Ni -•  LlNLtOUTFlLE); 


--  now  subtract  the  doM;  (overheat)  loop  aoiC  tcpot:  im  r«s..ts 

rUT(uLrTFI_E.  "FastPSN  executed  in  *5; 

FLOAT  10. FUT(OUTFIL£.  f loa t ( ( f as t pnn  elapsed 

dusay  elapseoi loe ) )  floa:(loop  count)); 

tex  t_io.  put_l  ine  (OlTTFllX,  "  eicv oscconds .  *  ) ; 

Pirr_LINE(  OUT  FILE,  "Returned  results"); 

printarray  (OUTFILE,  results,  dim,  nentrds); 

CLOSE (OUT FILE): 

end ; 

end  main; 
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with  Data_Struct_Pkg;  use  Data_Struct  Pkg; 
with  Build_Pkg; 
with  MergeDownPkg; 
with  Destroy_Pkg; 

procedure  FastPNN  (means  :  in  out  means_array_type; 

weights  :  in  out  weights_aLray_type, 

count  :  in  integer; 

dim  :  in  integer; 

ncntrds  :  in  integer; 

centroids  :  out  meansarraytype)  is 


Procedure  name:  FastPNN 

Purpose:  Outermost  routine  for  the  Pairwise  Nearest  Neighbor 
Clustering  algorithm 

Input  arguments: 

means  -  sample  point  array 

weights  -  sample  weight  array 

count  -  number  of  samples 

dim  -  dimensionality  of  the  sample  data 

ncntrds  -  target  number  of  clusters  to  form  the  data 

Output  arguments 

centroids  -  cluster  vector  results  from  FASTPNN  algorithm 


entryy  ;  kdent ry_ptr_type; 
tree  :  kdtree_ptr_type; 
k  :  integer; 

begin 

Build_Pkg. BuildKDtree(means,  weights,  count,  dim,  tree); 
HergeDown_Pkg - MergeDownKDt ree( tree , ncn  t  rds ) ; 

entryy  tree. root, bucket. Iists_array(0) ; 

k  0; 

for  i  in  0  . .  ncntrds  -  1  loop 
for  j  in  0  . .  dim  -  1  loop 

cencioids(k)  :=  entryy. mean_array( j ) ; 
k  :=•  k  +  1; 
end  loop; 

entryy  :»  entryy . next_array(0) ; 
end  loop; 

Des  t  r oy_Pkg . Fes  t  royKD  t  r ee ( tree ) ; 
end  FastPNN; 
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package  Data  Struct_Pkg  is 


—  This  package  specification  contains  the  declarations  of 
data  types  used  by  all  modules  in  the  FASTPNN  routine 


TRUEE  :  constant  ;=  -1; 

FALSEE  :  constant  :»  0; 

KDNODEE  :  constant  :*  0; 

KDBUCKETT  :  constant  :=  1; 

BUCKETSI2E  :  constant  8; 

KDMERGE  :  constant  :=  0.5; 

dim  :  constant  :=  2; 

—  Data  Structures  defined  in  the  main  program 

type  means  arraytype  is  array  (integer  range  <>)  of  float; 
type  veights_array_type  is  array  (integer  range  <>)  of  float; 

—  data  structure  used  in  CollapseKDnode 
type  integer_ptr_type  is  access  integer; 

—  data  structures  needed  in  indxx 

type  float_array_type  is  array  (integer  range  <>)  of  float; 
type  integer_array_type  is  array  (integer  range  <>)  of  integer; 
type  Kdentry; 

type  Kdentry_ptr_type  is  access  Kdentry; 

type  Kdentry_ptr_array_type  is  array  (integer  range  <>) 

of  Kdentry_ptr_type; 

type  meanarraytype  is  array  (integer  range  <>)  of  float; 
type  vmean_array_type  is  array  (integer  range  <>)  of  float; 
type  wsqran_array_type  is  array  (integer  range  <>)  of  float; 
type  wvar_array_type  is  array  (integer  range  <>)  of  float; 

type  Kdentry  is 
record 

next_array  :  Kdentry _ptr_arr3y_type(0  dim  -  1); 

splitleft  ;  integer; 

weight  :  float; 

mean_array  :  mean_array_type(0  ..  dim  -  1); 

vraeanarray  :  wmean_array_type(0  ..  dim  -  1); 

wsqmn  array  :  wsqmn  array  type(0  ..  dim  -  1); 

end  record; 

type  data_structure_type  is  (kdnode_type,  kdbuckettype) ; 
type  kdelera(data_structure  :  data_structure_type) ; 
type  kdelem_ptr_type  is  access  kdelem; 


—  Data  structure  kdelem_plr_arraytype  used  in  BalanceKDtree 
type  kdeiemptrar ray  type  is  array  (integer  range  <>) 
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of  kdelem_ptr_type 

type  kdnode  is 
record 

dindx  :  integer; 

lover  :  kdelem_ptr_type; 

upper  2  kdelem  ptrtype; 

end  record; 

type  kdbucket  is 
record 

count  :  integer; 

lists_array  :  Kdentry_ptr_array_type(0  ..  dim  -  1); 
entrya  :  Kdentry_ptr_type; 

entryb  ;  Kdentry_ptr_type; 

distort  ;  float; 

end  record; 

type  kdbucket_ptr_type  is  access  kdbucket; 

type  kdelem  (data_structure  :  data_structure_type)  is 
record 

typee  ;  integer; 

case  data_structure  is 
vhen  kdnode_type  »> 

node  2  kdnode; 

vhen  kdbucket_type  *> 

bucket  :  kdbucket; 

end  case; 
end  record; 

norb  n  2  kdelem(kdnode_type); 
norb  b  2  kdeleia(kdbucket_type) ; 

type  kdtree  is 
record 

dim  2  integer; 

root  2  kdelem_ptr_type; 

nbuckets  2  integer; 
nentries  2  integer; 
end  record; 

type  kdtree_ptr_type  is  access  kdtree; 
end  Data  Struct  Pkg; 
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with  Data_Stfuct_Pkg;  use  Data_Struct_Pkg; 
package  Rui''d_Pkg  is 


—  Package  Specification  Build_Pkg  contains  the  declaration  of 

—  all  subprograms  which  are  available  for  building  and 

—  Initialization  of  the  K-d  data  structure. 


procedure  BuildKDtree  (means  :  in  out  means_array_type; 

weights  :  in  out  weights_array_type; 

count  :  in  integer; 

dim  ;  in  integer; 

tree  :  out  Kdtree_ptr_type> ; 

function  CreateKDbucket  (dim  ;  integer)  return  kdeiem_ptr_type; 
end  Build  Fkg; 
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package  body  Build_Pkg  is 

function  CreateKDbucket  (dim  :  integer)  return  kdelem_ptr_type 

is  separate; 

procedure  Indxx  (arrin  :  in  float_array_type; 

indx  :  out  integer_array_type; 
n  :  in  integer)  is  separate; 

procedure  SortBucket  (elew  :  in  kdelem_ptr_type; 

dim  :  in  integer)  is  separate; 

function  CreateKDentry  (dim  :  in  integer)  return  kdentry_ptr_type 

is  separate; 

procedure  CreateFirstBucket  (means  :  in  out  tneans_array_type; 

weights  :  in  out  weights_array_type; 
count  :  in  integer; 
dim  :  in  integer; 
bucket  :  out  kdelem_ptr_type) 

is  separ? 

function  CreateKDtree  (dim  ;  integer)  return  kdtreeptrtype 

is  separate; 

procedure  BuildKDtree  (means  :  in  out  means_array_type; 

weights  :  in  out  welihts_array_type; 
count  !  in  integ^^) 
dim  s  in  integer; 

tree  ;  out  Kdtree_ptr_type)  is  separate; 

end  Build_Pkg; 
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separate  (Build_Pkg) 

procedure  BuildKDtree  (means  :  In  out  means_array  type; 

weights  :  in  out  weights_array_type; 

count  ;  in  integer; 

dim  :  in  integer; 

tree  :  out  Kdtrce_ptr_type)  is 


Procedure  name:  Buildkdtree 

Purpose:  Constructs  an  initial  KD  tree  from  the  sample  data 


Input  arguments: 

means  -  sample  point  array 

weights  -  sample  weight  array 

count  -  number  of  samples 

dim  -  dimensionality  of  the  sample  data 

Output  arguments: 

pointer  to  a  kdtree 


temp  tree  :  kdtree_ptr_type; 
bucket  :  kdelem_ptr_type; 

begin  —  BuildKDtree 

t6mp_tree  :*  CreateKDtree  (dim); 

CreatePirstBucket  (means,  weights,  count,  dim,  bucket); 
temp_tree.root  :»  bucket; 
temp_tree.nbuckets  :«  1; 
temp_tree.nentrias  count; 
tree  temp_trec; 
end  BuildKDtree; 
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separate  (Build_Pltg) 

procedure  CreateFirstBucket  (means  :  in  out  means_array_type; 

weights  :  in  out  weights_array_type; 

count  :  in  integer; 

dim  :  in  integer; 

bucket  ;  out  kdelem_ptr  type)  is 


—  Procedure  naxe: 

CreateFirstBucket 


—  Purpose: 

Create  and  initialize  the  first  bucket  in  a  KD  tree 


—  Input  arguments: 

means  -  sample  point  array 

weights  -  sample  weight  array 

count  -  number  of  samples 

dim  -  dimensionality  of  the  sample  data 


—  Output  arguments; 

bucket  -  pointer  to  first  bucket  in  the  Kdtree 


mean,  wvar  ;  mean_array_type(0  ..  dim  -1); 
kdptr,  Iptr  ;  kdentry  ptr  type; 

KDPTRO,  KDPTRl  :  KUENT RY_PTR_TYPE ; 
k  ;  integer; 

tempbucket  ;  Kdelem_ptr_type; 

begin  —  CreateFirstBucket 

tempbucket  CreateKDbucket(dim); 
tempbucket. bucket. count  count; 

k  ;■  0; 

tor  1  in  0  ..  (count  -  1)  loop 
kdptr  ;■  CreateKDentry(dim) : 
kdptr, weight  weights(i); 
for  j  in  0  ..  (dim  -  1)  loop 

kdptr .mean_array(j )  ;«  means(k); 

kdptr.vmean_array( j)  :«  kdptr .mean_array(j )  *  kdptr. weight; 
kdptr. wsqmn_array(j)  kdptr. mean_array(j )  * 

kdptr .wroean_array(j ) ; 

k  :>  k  1; 
end  loop; 
if  (i«0)  then 

tempbucket. bucket. Iists_array(0)  :»  kdptr; 

Iptr  :«  kdptr; 

else 

lptr.next_array(0)  ;»  kdptr; 

Iptr  kdptr; 
end  if; 
end  loop; 

SortBucket  (tempbucket,  dim); 
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bucket  tempBucket; 
end  CreateFirstBucket ; 
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with  Data_Struct_Pkg;  use  Data_Struct_Fltg; 
separate  (Build_Pkg) 

function  CreateObucket  (dim  ;  Integer)  return  kdelem_ptr_type  is 


—  Function  name: 

CreateKDbucket 


—  Purpose: 

Create  a  bucket  for  a  KD  tree 

—  Input  arguments 

dim  -  dimensionality  of  the  sample  data 

—  Output  arguments 

None 


~  -  Returns: 

—  returns  pointer  to  new  created  bucket 


bucket  ;  kdelem  ptrtype; 
begin 

bucket  :■  new  kdelem(kdbucket_type); 
bucket . typee  : »  KDBUCKETT ; 
return  bucket} 

end  CreateKDbucket; 
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separate  (Build_Pkg) 

function  CreateKDentry  (dim  :  in  integer)  return  kdentry_ptr_type  is 


—  Function  name: 

CreateKDentry 

—  Purpose: 

Create  a  KD  tree  bucket  entry  for  holding  a  sample  point 

—  Input  arguments; 

dim  -  dimensionality  of  the  sample  data 

—  Output  arguments: 

None 

—  Returns:  ' 

pointer  to  a  bucket  entry 


entryy  ;  kdentry_ptr_type; 
begin 

entryy  :■  new  Kdentry; 
return  entryy; 

end  CreateKDentry; 
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separate  (Bulld_Fkg) 

function  CreateKDtree  (dim  :  integer)  return  ltdtree_ptr_type  is 


—  Function  name: 

CreateKDtree 

—  Purpose: 

dynamically  allocates  storage  space  for  a  KD  tree  data  structure 

—  Input  Arguments: 

dim  -  dimensionality  of  the  tree  data 

—  Output  arguments 

None 

—  Returns: 

pointer  to  a  KDtree 


tree  :  kdtree_ptr_type; 
begin 

tree  :»  new  kdtree; 
tree. dim  dim; 
return  tree; 

end  CreateKDtree; 
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separate  (Build_Pkg) 

procedure  Indxx  (arrin  :  in  f loat_array_typG; 

indx  :  out  integer_array_typG; 
n  :  in  integer)  is 


—  Procedure  name:  Indxx 

Purpose;  Sort  an  array  of  Indeces  indx  based  on  the  data  arrin 
using  an  indexed  version  of  a  heap  sort.  Modified  from 
Numerical  Recipes  indexx 

Input  Arguments: 

arrin  -  array  of  data  used  for  sorting 
n  -  number  of  entries  in  arrin 

Output  arguments: 

indx  -  Indices  specifiying  the  order  of  data  in  arrin 


j,  i,  1,  ir,  indxt  :  integer; 
q  :  float; 

indxtemp  :  integer_array_type(0. .n-1) ; 

indx  j ,  indx_j_plus_l  :  integer; 
arrin_j ,  arrin_j_plus_l  :  float; 


bAfr  i  r\ 

for  k  in  0  . .  n-1  loop 
indxtemp(k)  k; 
end  loop; 

1  :«  n/2; 
ir  :-  n  -  1; 

while  TRUE  loop 

if  1  >  0  then 
1  1-  1; 

indxt  ;=*  indxtemp(l); 
q  :»  arrin(indxt); 

else 

indxt  :»  indxtempf ir) ; 
q  :=  arrin(lndxt) ; 
indxtemp(ir)  indxtemp(O); 
ir  ;»  ir  -  1; 
if  ir  =  0  then 

indxtemp(O)  :»  indxt; 
exit;  —  exit 

end  it; 
end  if; 

i  1; 

j  :=  2  *  (1  +  1)  -  1; 
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while  (j  <»  ir)  loop 

if  (j  <  ir  and  then  arrin( indxteiiip( j ) )  < 

arrin(indxteinp(  j  + 

i  i  *  li 
end  ifj 

if  (q  <  arrin(indxt€mp(j)))  then 
indxtenip(i)  indxtempC  j )  ; 

1  ii 

j  j  ^  1  +  1; 

else 

j  ir  +  1; 

end  if; 
end  loop; 

indxtemp(i)  ;•=  indxt; 
end  loop; 
indx  indxtemp; 

end  Indxx; 


))  then 
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separate  (Build_Pkg) 

procedure  SortBucket  (elem  :  in  kdelein  ptr_type; 

dim  :  in  integer)  is 


—  Procedure  name:  Sortbucket 

—  Purpose:  Sort  the  entries  in  a  KD  tree  bucket  across 

each  dimension  separately 

—  Input  arguments: 

elem  -  pointer  to  the  bucket  containing  the  entries 
to  be  sorted 

dim  -  dimensionality  of  the  sample  data 

—  Output  arguments:  None 


entryy  :  kdentry_ptr_type; 

epbffr  :  kdentiy_ptr_array_type(0  ..  elem. bucket. count-1) ; 
mnbffr  :  float  array_type(0  ..  elem. bucket . count-1 ) ; 
inbffr  :  integer_array_type(0  ..  elem. bucket .count-1) ; 
count  :  integer; 

begin 

count  elem. bucket. count; 
for  i  in  0  . ,  dim  -  1  loop 

entryy  ;«  elem. bucket. Iists_array(0): 
for  j  in  0  . .  count-1  loop 

ranbffr(j)  :*  entryy. mean_artay(i) ; 
epbffr(j)  entryy;  ~ 
entryy  :-  entryy .next_ar ray (0) ; 
end  loop; 

Indxxfmnbffr,  inbffr,  count); 

elem.bucket . lists_array(i)  :»  epbffr(inbffr(0)); 
for  j  in  1  . .  count  -  1  loop 

epbf f r( inbf f r( j  -  1) ) .next_array(l) 

epbf frvinbf f r( j )) ; 

end  loop; 

epbf ft(inbf fr(count  -  1) ) .next_array(i)  :»  null; 
end  loop; 

end  SortBucket; 
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with  Data_Struct_Pkg;  use  Data_Struct_Pkg; 
package  MergeOovn_Pkg  is 


—  Package  specification  6uild_Pkg  contains  the  declarations 

—  of  all  subprograms  which  are  available  for  reducing  the  K-d 

—  tree  built  in  module  BuildKDtree  to  the  the  specified  number 

—  of  output  centroids 


procedure  HergeDovnKDtree  (tree  :  in  kdtree_ptr_type; 

ncntrds  :  in  integer); 

end  MergeDown_pkg ; 
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with  Destroy_Pkg;  use  Destroy_Pkgj 
with  Build_Pkg;  use  Build_Pkg; 
package  body  MergeDownPkg  is 


procedure  ColiapseKDnode  (elem  ;  in  out  kdelem_ptr_type; 

dim  :  in  integer; 

bctr  :  in  out  integer)  is  separate; 

procedure  GetBucketStats  (elem  :  in  kdelem_ptr_type; 

dim  :  in  integer; 

mean  :  in  out  mean_array_type; 

wvar  :  in  out  wvar_array_type)  is 

separate; 

function  CreateKDnode  return  kdelem_ptr_type  is  separate; 

procedure  SplitBucket  (oldbucket  :  in  out  kdelem_ptr_type; 

dim  :  in  integer; 
bctr  ;  in  out  integer; 

bcktarray  ;  in  out  kdelem_ptr_array_type; 

bptr  :  in  out  integer; 

mean  :  in  out  mean_array_type; 

wvar  ;  in  out  wvar~array_type)  is 

separate; 

procedure  BalanceKDtree  (tree  ;  in  kdtree_ptr_type; 

bcktarray  :  in  out 

kdelem_ptr_array_type)  is  separate; 

procedure  AssessCandidate  (elem  :  in  kdelem_ptr_type; 

dim  :  in  integer; 

barray  :  in  out  kdelera_ptr_array_type; 

ncount  :  in  out  integer)  is  separate; 

function  BucketCompare  (one  :  Kdelem_ptr_type; 

two  ;  Kdelem_ptr_type) 

return  boolean  is  separate; 

procedure  Quicksort  (bcktarray  ;  in  out  Kdelem_ptr_array_type) 

is  separate; 


procedure  ReduceKDbucket  (elem  :  in  kdelemptrtype; 

dim  :  in  integer)  is  separate; 

procedure  CompressKDtree  (tree  :  in  kdtreeptrtype; 

nmerge  :  in  integer; 
bcktarray  ;  in  out 
kdelem_ptr_array_type)  is  separate; 

procedure  MergeDownKDtree  (tree  ;  in  kdtreeptrtype; 

ncntrds  :  in  integer)  is  separate; 

end  MergeDown_Pkg; 
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separate  (MergeDovn_Pkg) 

procedure  AssessCandidate  (elem  :  In  kdelein_ptr_type;  dim  :  in  integer; 

barray  :  in  out  k,delem_ptr_airay_type; 
ncount  :  in  out  integer)  is 


—  Procedure  name:  Assesscandidate 

—  Purpose:  Determine  the  pair  u£  bucket  entries  which  produces 

the  minimal  mean  distortion  when  merged 

—  Input  arguments: 

elem  -  pointer  to  the  KD  tree  bucket  under  evaluation 
barray  -  array  of  pointer  to  all  buckets  in  the  KD  tree 
ncount  -  counter  initialized  to  zero 

—  Output  arguments: 

barray  -■  array  of  pointers  to  all  buckets  that  can  be  merged 
ncount  -  counter  used  to  keep  track  of  the  the  total  number 
of  buckets  that  can  be  merged 


ientry,  j entry  :  kdentry_ptr_type; 
reduction,  dotprd,  diff  :  float; 
firsttime  :  integer; 

begin 

first time  TRUEE; 

if  (elem. bucket. count  >1)  then 

ientry  :=  elem. bucket. Iists_array(0) ; 
for  i  in  0  ..  elem. bucket. count  -  1  loop 
jentry  :=  ientry. next_array(0) ; 
for  j  in  i  +  1  . .  elem. bucket. count  -  1  loop 
dotprd  :=  0.0; 
for  k  in  0  ..  dim  -  1  loop 

diff  ;=  ientry. roean_array(k)  -  jentry. mean_array(k) 
dotprd  :»  dotprd  diff  *  diff; 
end  loop; 

reduction  :=  (dotprd  *  ientry. weight  *  jentry. weight)  / 
(ientry. weight  +  jentry.weight) ; 

if  ((reduction  <  elem. bucket. distort)  or 
(firsttime  *  TRUEE))  then 
elem. bucket. distort  reduction; 

elem. bucket. entiya  ;=  ientry; 
elem. bucket. entryb  ;=  jentry; 
firsttime  FALSER; 
end  if; 

jentry  ;=  j entry, next_array(0) ; 
end  loop;  —  end  j  loop 
ientry  ;=  ientry. next_array(0) ; 
end  loop;  —  end  i  loop 

barray(ncount)  :=  elem; 
ncount  :=  ncount  +  1; 


B-21 


end  if 


end  AssessCandldate; 
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separate  (MergeDovn_Pk£) 

procedure  BalanceKDtcee  (tree  :  in  kdtree_ptr_type; 

bcktarray  ;  in  out  kdelein_ptr_array_type)  is 


—  Procedure  name:  BalanceKDtree 

—  Purpose:  Redistribute  the  entries  in  a  KD  tree  so  that  each  bucket 

has  approximately  the  same  n'lmber  of  entries 

—  Input  arguments: 

tree  -  pointer  to  a  KD  tree 

—  Output  arguments: 

bcktarray  -  array  used  to  retain  a  pointer  to  each  bucket  after 
balancing 


bptr  :  integer; 

mean  ;  mean_array_type(0  ..  dim  -  1)  :=  (others  *>  0.0); 
war  :  war_array_type(0  ..  dim  -  1)  :»  (others  *>  0.0); 

begin 

if  (tree.root.typee  «  KDNODEE)  then 

CollapseK0node( tree. root >  tree. dim,  tree.nbuckets); 
end  if; 

bcktarray(O)  :*  tree. root; 
bptr  :«  1; 

SplitBucket( tree. root,  tree. dim,  tree.nbuckets,  bcktarray, 
bptr,  mean,  war); 

end  BalanceKDtree; 
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separate  (MergeDown_Pkg) 

function  BucketCompare  (one  :  Kdelem_ptr_type;  two  :  kdelem  ptr  type) 

return  boolean  i 

begin 

if  one. bucket. distort  <  two. bucket. distort  then 
return  FALSE; 

else 

return  TRUE; 
end  if; 

end  BucketCompare; 
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separate  (MergeOovn_Pkg) 

procedure  CollapseKDnode  (elem  :  in  out  kdelem_ptr_^type; 

dim  :  in  integer; 

bctr  :  in  out  integer)  is 

FIRSTIME  :  BOOLEAN  TRUE; 

ptr,entptr,  lentry,  rentry,  testentry;  kdentry_ptr_type; 
bucket  :  kdelem_ptr_type; 

begin 

if  elem. node. lower. ty pee  /■  KDBUCKETT  then 

CollapseKDnode  (elem. node. lover ,  dim,  bctr); 
end  if; 

if  elem. node. upper. ty pee  /»  KDBUCKETT  then 

CollapseKDnode  (elem. node. upper ,  dim,  bctr); 
end  if; 

bucket  elem. node. lover ; 

for  1  in  0  . .  dim  -  1  loop 

entptr  :»  elem. node. lower. bucket. lists_array(i); 
lentry  ;■  elem. node. lover. bucket. lists_array(i); 
rentry  clem. node. upper. bucket. lists_array(i); 

testentry  i«  elem. node. upper. bucket. lists_array(l) ; 

while  (lentry  /»  null)  and  (rentry  /■  null)  loop 

if  lentry. mean_array(i)  <  rentry.mean_array(i)  then 

if  (FIRSTIME  -  TRUE)  then 

elem. node. lower. bucket. lists_array(i)  lentry; 

entptr  :=  lentry; 

FIRSTIME  FALSE; 
else  —  FIRSTIME  -  FALSE 

entptr. next_array(i)  ;»  lentry; 
entptr  :«  lentry; 
end  if; 

lentry  ;»  lentry.next_array(i); 

else  —  case  where  leRtry.mean_array(i)  >- 

rentry. sean_ar ray ( i ) 


if  (FIRSTIME  -  TRUE)  then 

elem. node. lover .bucket . lists  array(i)  ;»  rentry; 
entptr  ;=  rentry;  ~ 

FIRSTIME  :«=  FALSE; 

else 

entptr .next_array(i)  :=  rentry; 
entptr  :=.  rentry; 
end  if; 

rentry  :*  rentry.next_array(i) ; 
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end  loop; 


if  (lentry  /«  null)  then 

while  (lentry  /.  null)  loop 

entptr.next_array(i)  :•  lentry; 
entptr  lentry; 
lentry  :=  lentry.next_array(i) ; 
end  loop; 
end  if; 

if  (rentry  /»  null)  then 

while  (rentry  /■  null)  loop 

entptr. next_array(i)  :»  rentry; 
entptr  :«  rentry; 
rentry  ;■  rentry.next_array(i) ; 
end  loop; 
end  if; 

FIRSTIME  TRUE; 
end  loop; 


bucket. bucket. count  :=  elem. node. lover. bucket. count  + 

elem. node . upper . bucke  t . coun  t 

DestroyRObucket  (elem. node. upper) ; 

DestroyKDnode  (elem); 

bctr  :»  bctr  -  1; 
eles  bucket; 

end  CollapseKDnode; 
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SI 


separate  (MergeDown_Pkg) 

procedure  CompressKDtree  (tree  :  in  kdtree_ptr_type;  nmerge  :  in  integer; 

bcktarray  :  in  out  kdelem_ptr  array_type)  is 


—  Procedure  name;  CompressKDtree 

—  Purpose:  Proceduie  used  to  merge  bucket  entry  pairs  into  single  bucket 

entries  for  a  fixed  fraction  of  the  total  number  of  buckets 

—  Input  arguments; 

tree  -  pointer  to  a  KD  tree  undergoing  the  merge 

nmerge  -  number  of  bucket  pairs  to  merge 

bcktarray  -  array  of  pointers  to  all  buckets  in  the  tree 

—  Output  arguments; 

none 


ptr  ;  kdentry_ptr_type; 
nco'int  ;  integer; 

TESTPTR  ;  kdentry_pt retype; 
nmerge  temp  ;  integer  ;-  nmerge; 

begin 

ncount  ;■  0; 

—  for  each  bucket,  one  at  a  time,  find  the  candidate  pair, 

—  ENTRYa  and  ENTRYB,  which  has  the  least  weighted  square 

—  error  (DISTORT)  —  between  them 

for  i  in  0  ..  tree.nbuckets  -  1  loop 

AssessCandidate(bcktarray<i),  tree. dim,  bcktarray,  ncount); 
end  loop; 

if  (nmerge  >  ncount)  then 
nmerge  temp  ncount; 
end  if; 

if  (ncount  >  1)  then 

Quicksort(bcktarray(0  ..  ncount-l)); 
end  if; 

for  i  in  0  . .  nmerge  temp  -  1  loop 

ReduceKDbucket(bcktarray(i) ,  tree. dim) ; 
end  loop; 

tree.ncntries  ;«  tree.nentries  -  nmergetemp; 
end  CompressKDtree; 
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separate  (HergeDovn_Pkg) 

function  CreateKDnode  return  k.delein_ptr_type  is 


—  Function  name:  CreateKDnode 

—  Purpose:  Create  a  KD  tree  node 

—  Input  arguments: 

none 

—  Returns: 

pointer  to  the  newly  created  KD  tree  node 


node  :  kdelem_ptr_type; 
begin 

node  :=  new  kdelera(kdnode_typs) ; 
node-typee  ;»  KDNODEE; 
return  node; 

end  CreateKDnode; 


separate  (MergedownPkg) 

procedure  GetBucketStats  (elem  ;  in  kdelem_ptr_type; 

dim  :  in  integer; 

mean  :  in  out  inean_array_type ; 

war  ;  in  out  wvar_array_type)  is 

entryy  ;  kdentry  ptrtype; 
wgtsum  :  float; 

begin 

vgtsum  0.0; 

for  i  in  0  . .  dim  -  1  loop 
mean(l)  0.0; 
vfvar(i)  :»  0.0; 
end  loop; 

entryy  ;■  elem. bucket . Iists_array(0) ; 
while  (entryy  /»  null)  loop 

vgtsum  :=  wgtsum  +  entryy. weight; 
for  i  in  0  . .  dim  -  1  loop 

mean(l)  :»  mean(i)  +  entryy.wmean_array(i); 
wvar(i)  ;»  wvar(i)  +  entryy. wsqmn_array(i) ; 
end  loop; 

entryy  ;•  entryy. next_array(0); 
end  loop; 

for  i  in  0  . .  dim  -1  loop 

mean(i)  me8n(i)  /  wgtsum; 

war(i)  :■  (wvar(i)  /  wgtsum)  -  (mean(i)  *  mean(l)); 
end  loop; 

end  GetBucketStats; 
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separate  (Mergedovn_Pkg) 

procedure  MergeDovnKDtree  (tree  :  In  k<ltree_ptr_type; 

ncntrds  :  in  integer)  is 


•—  Procedure  name;  MergeDovnKDtree 

—  Purpose:  Reduce  a  KD  tree  to  a  single  bucket  having  ncntrds 

entries  using  the  PNN  algorithm 

—  Input  arguments: 

tree  -  pointer  to  a  KD  tree 

ncntrds  -  desired  number  of  entries  after  merging 

—  Output  arguments; 

None 

—  Returns: 

Nothing 


ptr  :  kdentry_ptr_type; 
nmerge,  ntile  :  integer; 

maxbkts  ;  integer  :■  tree.nentries  /(BUCKETSIZE  /  2); 
bcktarray  :  kdelem_ptr_array_type(0  ..  maxbkts); 

begin 

while  (tree.nentries  >  ncntrds)  loop 

BalanceKDtree  (tree,  bcktarray); 

if  KDHERGE  *  float ( tree. nbuckets)  >  1.0  then 

—  truncate  by  the  integer  conversion  by  subtracting  .5 
ntile  :»  integer ( (KDHERGE  *  float( tree. nbuckets))  -  0.5) 


else 

ntile  :«  1; 
end  if; 

if  (tree.nentries  -  ncntrds)  <  ntile  tnen 
nmerge  tree.nentries  -  ncntrds; 

else 

nmerge  ntile; 
end  if; 

CompressKDtree  (tree,  nmerge,  bcktarray); 
end  loop; 

if  tree.root. typee  «  KDNODEE  then 

CollapseKDnode  (tree.root,  tree. dim,  tree.nbuckets) ; 
end  if; 

end  MergeDovnKDtree; 
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separate  (MergedownPkg) 

procedure  Quicksort  (bcktarray  ;  in  out  kdeleiD_ptr_array_type)  is 

SIZE  :  INTEGER  bcktarray' length; 

FIRST  :  INTEGER  bcktarray ' FIRST; 

LAST  :  INTEGER  bcktarray' LAST; 

P  :  INTEGER; 

procedure  SWAP(A,  B  :  in  out  KDELEM_PTR_TYPE )  is 
TEMP  :  KDELEM_PTR_TYPE ; 
begin 

TEMP  A; 

A  :  •  B ; 

B  TEMP; 
end  SWAP; 

procedure  PARTITION{ bcktarray  :  in  out  KDELEM_PTR_ARRAY  TYPE; 

NEWPIVOT  :  out  INTEGER)  is 

PIVOTVALUE  :  KDELEM_PTR_TYFE ; 

FIRST  ;  INTEGER  ;=  bcktarray' FIRST; 

LAST  :  INTEGER  bcktarray' LAST; 

LOVER  :  INTEGER; 

UPPER  :  IirrEGER; 

LASTOPEN  :  INTEGER; 

begin 

LOVER  ;=  FIRST; 

UPPER  LAST; 

PIVOTVALUE  bcktarray(FIRST) ; 

USTOPEN  :=  FIRST; 

while  LOVER  <  UPPER  loop 

while  UPPER  >  LOVER  and  then 

BucketCompare(bcktarray(upper) ,  PIVOTVALUE)  loop 

UPPER  UPPER  -  1; 

end  loop; 

bcktarray (LASTOPEN)  ;»  bcktarray(UPPER) ; 

LASTOPEN  :=  UPPER; 

while  LOVER  <  UPPER  and  then 

Bucke tCompare( PI VOTVALUE , bcktarray ( lower ) )  loop 

LOVER  LOVER  +  1; 

end  loop; 

bcktarray(LASTOPEN)  :=•  bcktarray( LOVER) ; 

LASTOPEN  :=  LOVER; 
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end  loop; 


bcktarray(LASTOPEN)  !-  PIVOTVALUE; 
NEWPIVOT  LASTOPEN; 
end  PARTITION; 

begin  QUICKSORT 

if  SIZE  <=  1  then 
null; 
else 

PARTITICN(bcktarray , P) ; 
QUICKSORT(bcktarray( FIRST  ..  P  -  1)); 
QUIC^’SORT(bcktarray(P  +  1  ..  LAST)); 
end  if; 

end  Quicksort; 
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with  Data_Struct_Pk;g;  use  Data_Struct_Pkg; 

separate  '(MergeDown_Pkg)  ~ 

procedure  ReduceKObucket  (elem  :  in  kdelem_ptr_type; 

dim  :  in  integer)  is 


—  Procedure  name:  ReduceKDbucket 

—  Purpose:  Merges  a  pair  of  bucket  entries  into  a  single  entry 

—  Input  arguments: 

elem  -  pointer  to  the  KO  tree  bucket  whose  entries  are  to 
be  merged 

dim  ~  dimension  of  the  data  within  the  entries 

—  Output  arguments: 

None 


last_entry,  ientry,  jentry,  leptr,  oldptr  :  kdentry_ptr_type; 
newweight  :  float; 
rmcnt  :  integer; 

R£MOVE_FIRST  ENTRY  :  BOOLEAN  :=.  TRUt.; 

FIRST_ENTRY  7  BOOLEAN  TRUE; 
ptr,  TEST_PTR  :  kdentry  ptr_type; 

begin 

ientry  :>  elem. bucket. entrya; 
jentry  elem.bucket.entryb; 

—  remove  ientry  and  jentry  from  the  list 
for  i  in  0  . .  dim  -  1  loop 

leptr  :=  elem. bucket. llsts_array(i) ; 

rmcnt  ;«>  0; 

while  rmcnt  <  2  loop 

if  ((leptr  /»  ientry)  and  (leptr  /«  jentry))  then 
last_entry  leptr; 
leptr  :*>  leptr.next_array(i) ; 

REMOVE_FIRST_ENTRY  :-  FALSE; 

else 

if  REMOVE_FIRST_ENTRY  -  TRUE  then 

elem. bucket. lists_array(i)  :«  leptr.next_array(i) ; 
else 

last_entry.next_array(i)  :-  leptr. next_array( i) ; 
end  if; 

leptr  :=  leptr. next_array(i ) ; 
rmcnt  :=  rmcnt  +  1; 
end  if; 
end  loop; 

REMOVE_FIRST_ENTRY  :=  TRUE; 
end  loop; 

newweight  :=  ientry. weight  +  jentry. weight; 
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for  i  in  0  . •  dim  -  1  loop 

ientry.mean_array(i)  :=*  <lentry.mean_array(i)  * 

”  lentry. weight  +  jentry.ir.ean_array(i) 

*  jentry. weight)  /  newweight; 


ientry . wmean_ar ray ( i ) 
i en  t  ry . wsqmn_ar r ay ( i ) 
end  loop; 


:■  lentry. mean_array(i)  *  newweight; 

:=  ientry. mean_array(i)  * 

lentry. wmean  array(i); 


ientry  .weight  newweight; 

—  Reinsert  ientry  into  the  list  in  the  proper  order 


for  i  in  0  . .  dim  -  1  loop 

leptr  :=  elem. bucket. lists_array(i); 

—  traverse  entries  that  will  remain  unchanged  in 

—  original  list 

while  (leptr  /=  null)  and  then 

(leptr.mean_array(i)  <  ientry. mean_array(i))  loop 

oldptr  :*  leptr; 

leptr  ;=  leptr.next_array(i); 

FIRST_ENTRY  ;=  FALSE; 
end  loop; 

—  insert  ientry 

if  (FIRST_ENTRY  =  TRUE)  then 

elem. bucket. lists_array(i)  ientry; 
else 

oldptr.next_array(i)  :=  lentry; 
end  if; 

lentry. next_array(i)  :=  leptr; 

FIRST_ENTRY  TRUE; 
end  loop; 

elem. bucket. count  j«  elem. bucket. count  -  1; 


Des  t  r oyKDen  try (jentry); 
end  ReduceKDbucke t ; 
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with  Build_Pkg;  use  Build_Pkg;  —  to  "see"  CreateKDbucket 
separate  (HergeDown_Pkg) 

procedure  SplitBucket  (oldbucket  :  in  out  kdelem_ptr_type; 

diia  :  in  integer; 
bctr  ;  in  out  integer; 

bcktarray  :  in  out  kdelem_ptr_array_type; 

bptr  :  in  out  integer; 

mean  :  in  out  mean_array_type; 

war  :  in  out  wvar  array  type)  is 


—  Procedure  name:  Splitbucket 

—  Purpose:  Recursive  procedure  used  to  split  a  bucket  into  two 

smaller  buckets  having  half  as  many  entries 


j,  brount,  medindx  :  integer; 
newnode  :  kdelem_ptr_type; 
newbucket  :  kdelem_ptr_type; 
dummy  :  kdelem_ptr_type; 
oldptr,  newptr  :  kdentry_ptr_type; 
entryy  :  kdentty_ptr  type; 

FIRSTIME_LEFT  :  BOOLEAN  :=  TRUE; 

FIRSTIME_RIGHT  :  BOOLEAN  :=  TRUE; 

begin 

if  (oldbucket. bucket. count  >  BUCKETSIZE)  then 

GetBucketStats  (oldbucket,  dim,  mean,  war); 

—  find  dimenson  with  largest  variance 

j  0; 

for  i  in  1  . .  dim  -  1  loop 

if  (wvar(i)  >  wvar(j))  then 

j  :=  i;  —  j  is  dimenson  with  largest  variance 

end  if; 
end  loop; 

bcount  :«  oldbucket. bucket. count;  —  bcount  is  BUCKETCOUNT 
medindx  :a  (bcount  1)  /  2;  —  median  index 

newnode  :=  CreateKOnode; 

newnode. node. dindx  ;=  j;  —  dimension  to  be  split 
newbucket  ;=  CreateKDbucket (dim) ; 
newnode. node. lower  :=  oldbucket; 
newnode. node. upper  :=  newbucket; 

bctr  :=  bctr  +1;  —  increment  number  of  buckets  in  tree 


B-35 


bcktarray(bptr)  :=  ncvbucket; 

bptr  !»  bptr  +1;  —  increment  array  index  of  bcktairay 

traverse  the  entries  below  the  median  value  in  the 

—  j  dimension  (the  dimension  with  the  largest  variance) 

—  and  set  entryy.splltlef t  =  TRUE 

entryy  :=  oldbucket. bucket. lists_array(j); 
for  i  in  0  . .  medindx  -  1  loop 
entryy. splitleft  :*  TRUEE; 
entryy  :=  entryy.next_array(j ) ; 
end  loop; 

traverse  the  entries  above  the  median  value  in  the 

—  j  dimension  (the  dimension  with  the  largest  variance) 

—  and  set  entryy .splitleft  =  FALSE 

for  i  in  medindx  . .  bcount  -  1  loop 
entryy. splitleft  :=  FALSER; 
entryy  :=  entryy. next_array( j ) ; 
end  loop; 

oldbucket. bucket. count  ;=  medindx; 

newbucket. bucket .count  ;=  bcount  -  medindx; 

for  i  in  0  . .  dim  -  1  loop 

oldptr  j!»  oldbucket. bucket. lists_atray(i) ; 
newptr  ;■  newbucket. bucket. lists_array(i); 
entryy  :«  oldbucket. bucket. lists_array(i); 

while  (entryy  /=  null)  loop 

if  (entryy. splitleft  =  TRUEE  and 

FIRSTIME_LEFT  =  TRUE)  then 
oldbucket. bucket. lists_array(i)  entryy; 
oldptr  ;«  entryy; 

FIRSTIME_LEFT  FALSE; 

elsif  (entryy.splltlef t  =  TRUEE  and 

FIRSTIME_LEFT  -  FALSE)  then 
oldptr. next_array(i)  ;■  entryy; 
oldptr  entryy; 

elsif  (entryy. splitleft=FALSEE  and 

FIRSTIME_RIGHT=TRUE)  then 
newbucket. bucket. lists_array(i)  entryy; 
newptr  ;=  entryy; 

FIRSTIME_RIGHT  ;=  FALSE; 

elsif  (entryy.splltlef ts^FALSEE  and 

FIRSTIMERIGHT-FALSE)  then 
newptr.next_array(i)  entryy; 
newptr  ;=•  entryy; 
end  if; 
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entryy  enttyy.next_array(i); 
end  loop; 

FIRSTIME_LEFT  TRUE; 

FIRSTIME_RIGHT  TRUE; 

oldptr.next_array(i)  ;»  null; 
newptr.next_array(i)  :=  null; 
end  loop; 


SplitBucket(nevmude. node. lover,  dim,  bctr, 
bcktarray,  bptr,  mean,  war); 

SplltBucket(nevnode. node. upper,  dim,  bctr, 
bcktarray,  bptr,  mean,  war); 

oldbucket  ;=  newnode; 

end  if; 

end  SplitBucket; 
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with  Data_Struct_Pk*;  use  Data_Struct_Pkg; 
with  UNCHECKEDDEALLOCATION ; 
package  Destroy_Pkg  is 


—  This  package  specification  contains  the  dynamic  deallocation 
—  routines  which  are  visible  are  to  be  visible 

procedure  DestroyKDnode  (node  :  in  out  Kdelem_ptr_type) ; 

procedure  DestroyKDtree  (tree  !  in  out  Kdtree_ptr_type) ; 

procedure  DestroyKDbucket  (bucket  :  in  out  kdelem_ptr_type) ; 

procedure  DestroyKDentry  (entryy  :  in  out  kdentry_ptr_type) ; 

end  DestroyPkg; 
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package  body  Destroy_Pkg  is 


—  This  package  body  contains  the  routines  for  deallocated 

—  memory  which  was  previously  dynamically  allocated. 


procedure  DestroyKDnode  (node  ;  in  out  kdelem_ptr_type)  is 

separate; 

procedure  DestroyKDentry  (entryy  ;  in  out  kdentry_ptr_type)  is 

separate; 

procedure  DestroyKDbucket  (bucket  ;  in  out  kdelem_pt retype)  is 

separate; 

procedure  DestroyLastBucket  (bucket  :  in  out  kdelem_ptr_type)  is 

separate; 

procedure  DestroyKDtree  (tree  :  in  out  kdtL-ee_ptr_type)  is 

separate; 

end  Destroy_Pkg; 
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separate  (Destroy_Pkg) 

procedure  DestroyKDbucket  (bucket  :  in  out  Kdelem  ptr  type)  is 


—  Procedure  name:  DestroyKDbucket 

—  Purpose:  Destroy  a  Kd  tree  bucket 

—  Input  arguments: 

bucket  -  bucket  to  be  destroyed 

—  Output  arguments: 

bucket  -  bucket  is  null  upon  successful  execution 

of  FREE (bucket) 


procedure  FPvEE  is  new  UNCHECKED_DEALLOCATION(kdelem, 

kdelem  ptr_type); 


begin 

FREE (bucket); 
end  DestroyKDbucket; 
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separate  (Destroy _Pkg) 

procedure  DestroyKDentry  (entryy  ;  in  out  Kdentry_ptr_type)  is 


—  Procedure  name;  DestroyKDentry 
Purpose:  Destroy  a  KD  tree  bucket  entryy 
Input  arguments: 

entryy  -  pointer  to  the  entryy  to  be  destroyed 

—  Output  arguments; 

entryy  -  points  to  null  upon  successful  execution  of 

FREE(entryy) 


procedure  FREE  is  nev  UNCHECKED_DEALLOCATION  (Kdentry, 

Kdentry_ptr_type) ; 


begin 

FREE (entryy) ; 
end  DestroyKDentry; 


1. 


separate  (Destroy_Pkg) 

procedure  DestroyLastBucket  (bucket  ;  in  out  Kdelem_ptr_type)  is 


Procedure  name  :  DestroyLastBucket 

—  Purpose  :  Destroy  the  last  bucket  in  a  tree 

Input  arguments: 

bucket  -  bucket  to  be  destroyed 

—  Output  arguments: 

None 


entryy,  next  ;  kdentry_ptr_type; 
begin 

entryy  bucket. bucket. Iists_ariay(0) ; 
while  entryy  /-  null  loop 

next  ;«  entryy. next_array(0) ; 
DestroyKDentry  (entryy); 
entryy  :«  next; 
end  loop; 

Des  t  rovKDbucke  t ( bucke  t ) ; 
end  DestroyLastBucket; 
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separate  (Destroy^Pkg) 

procedure  DestroyKDnode  (node  :  in  out  Kdelem  ptr  type)  Is 


—  Procedure  name:  DestroyKDnode 

—  Purpose: 

Deallocate  a  ICO  tree  node 


—  Input  argum?.nts; 

node  -  pointer  to  the  entry  to  be  destroyed 

—  Output  arguments: 

node  -  point  to  null  upon  successful  execution  of  FREE(node) 


procedure  FREE  is  new  UNCHEClCED_DEALLOCATION(kdelem, 

kdelem  ptr  type); 


begin 

FREE(node' ; 
end  Df  s  t 'oyKDnode ; 
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separate  (Destroy_Pkg) 

procedure  DestroyKDtree  (tree  :  in  out  kdtree_ptr_type)  is 


—  Procedure  name:  DestroyKDtree 

—  Purpose:  Deallocate  a  KD  tree  structure 

Input  arguments; 

tree  -  pointer  to  the  KD  tree 

—  Output  arguments: 

tree  -  pointer  to  null  upon  successful  completion  of 
FREE (tree) 


procedure  FREE  is  new  UNCHECKED_DEALLOCATION(kdtree, 

kdtree  ptr  type); 


begin 

if  tree. root  /»  null  then 

DestroyLastBucket(tree.root) ; 
end  if; 

FREE(tree); 

end  DestroyKDtree; 
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package  timer  is 


—  See  bottom  for  comments  on  hov  to  use  the  timer 


type  microsec_timer  is  private) 
procedure  init_timer(  obj  :  out  microsec_limer  ); 
function  elapsed_time(  obj  ;  in  micro5ec_timer  )  return  integer; 
function  identity(  arg  :  integer  )  return  integer; 
function  alvays_true  return  boolean; 
private 

type  microsec_timer  is  new  integer; 
end  timer; 

—  Calling  sequence:  Call  the  subprograms  in  this  order: 

loop_count  :  constant  integer  100;  —  say... 

dummy_timer  :  timer. microsec  timer; 
dummy_elapsed_tiffle  ;  integerT 
dummy_arg  :  integer; 

— '  iiiy  tiwei  :  tinier  .microsec  timer; 

my”elapsed_tirae  :  integerj 

begin 


now  loop  through  dummy  loop  to  determine  the  length 
of  time  it  takes  to  run  through  the  loop  itself. 

Ue  will  subtract  this  out  later. 

timer . init_timer(  dummy^timer  ); 

for  i  in  1. . loop_count 

loop 

dummyarg  :«  timer. identity(  dununy_arg  ); 
end  loop; 

if  timer. always_true 
then 

dummy_elapsed_time  :*  timer. elapsed_time(  duramytiraer  ); 
end  if; 

—  now  run  through  the  sane  loop  and  add  the  actual  code  that 

—  you  want  to  time. 

timer . init_timer(  my_tiraer  ); 
for  i  in  1. . loopcount 
loop 

—  do  something  here  that  needs  timing. . . 
timer . identity(  dummyarg  ); 
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end  loop; 

if  timer .always_t rue 
then 

my_elapsed_tiroe  :»  timer. elapsed_time(  my_tiraer  ); 
end  if; 

—  now  subtract  off  the  dummy  loop  time  and  divide  by  the 

—  number  of  times  that  the  loop  occurs.  Thi.s  equals  one 

—  iteration  of  the  interesting  code. 

time_for_one_iteration 

(  my_elapsed_tirae  -  dummy_elapsed_time  )  /  loop_count  ); 

note  that  the  value  of  mytimer  is  not  changed  by  a  call  to 
timer .elapsed_time.  You  can  have  as  many  timers  as  you  wish. 
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with  system  ;  use  system  ; 

vith  condition  handling  ;  use  condition_handling  ; 
with  starlet  ;”use  starlet  ; 

package  body  timer  is 


—  Package:  timer 

>  —  Description:  This  package  provides  two  subprograms  which  are  called 

—  to  time  operations  in  the  vms  environment > 

—  References:  Please  see  the  VAX/VMS  System  Services  Reference  Manual 

—  under  the  $GETJPI  system  service  description. 


—  Revision  history:  S.  French  25'Jan-1989 


—  Host  processor:  VAX/VMS  4.7 


—  Package  dependencies:  This  package  interfaces  with  the  VMS  system 

—  service  $GETJPI  to  get  the  cpu  time. 


function  cpu_time^clock  return  integer  is 
cputim  :  integer  ; 
pragma  volatile  (  cputim  )  ; 
jpi_status  :  cond_value_type  ; 
jpi3item_list  :  constant  item_list_type  :■ 

~  ~  (  (  -^  »  jpi_cputim'",  cputim'address  ,  address_2ero  )  , 

(0,0  , "address  zero  ,  address  zero  )  )  ;” 


begin 

—  call  get j pi  to  set  cputim  to  total  accumulated  cpu  time 

—  (in  millisecond  tics) 

getjpi  (  status  ->  jpi_status  ,  Itmlst  ->  jpi_item_list  )  ; 
return  (cputim  *  10) j 
end  cpu_time_clock  ; 

procedure  init_timer(  obj  :  out  micro5ec_timer  )  is 
begin 

obj  :»  microsec_timer(  cpu_time_clock  ); 

*  end  inlt_timer;  " 

function  elapsed_time(  obj  :  in  microsec_timer  )  return  integer  is 
begin 

return  cpu_time_clock  -  integer(  obj  ); 
end  elap&ed_time; 

function  identity(  arg  :  integer  )  return  integer  is 
some  value  :  integer  0; 
begin 

some_value  ;«  some_value  +  arg; 
return  some  value; 
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end  identity 


function  alvays_true  return  boolean  is 
bool_value  :  boolean  true; 
begin 

return  bool_value; 
end  always_true; 

begin 

null; 

end  timer; 
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package  timer  is 


—  See  bottom  for  comments  on  hov  to  use  the  timer 


type  microsec_timer  is  private; 
procedure  init_tiner(  obj  :  out  microsec  timer  ); 
i  function  elapsedtime(  obj  :  in  microsec_timer  )  return  Integer; 

function  identity<  arg  :  integer  )  return  integer; 

i 

function  alvays_true  return  boolean; 
private 

type  raicrosec_timer  is  new  integer; 
end  timer; 


—  Package:  timer 

—  Description:  This  package  provides  two  subprograms  which  are  called 

—  to  time  operations. 

—  Calling  sequence:  Call  the  subprograms  in  this  order: 

loop_count  :  constant  integer  :«  100;  —  say... 

duramy_timer  :  timer. raicrosec_timer; 

dummy_elapsed_time  :  integer; 

duramy_arg  ;  integer; 

ray_tiraer  :  timer. microsec_timer; 

my_elapsed_time  :  integer; 

—  begin 

now  loop  through  dummy  loop  to  determine  the  length  of  time 
it  takes  to  run  through  the  loop  itself.  We  will  subtract 
this  out  later. 

timer, init_timer(  dummy_timer  ); 
for  i  in  1. .loop_count  loop 

dumrayarg  :»  timer. identity(  dummy _arg  ); 
end  loop; 

if  timer.alwaystrue  then 

dummy _elapsed_time  :=  timer. elapsed  time(  dummy  timer  ); 
end  if; 

now  run  through  the  same  loop  and  add  the  actual  code  that 
you  want  to  time. 

timer,  ini r._timer(  my  timer  ); 
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for  i  in  1. .loop_count  loop 

do  something  here  that  needs  timing. . . 
timer. identityC  dummy_arg  ); 
end  loop; 

if  timer. always_true  then 

my_elapsed_time  :=  timer. elapsedtime(  my_timer  ); 
end  if; 

now  subtract  off  the  dummy  loop  time  and  divide  by  the 
number  of  times  that  the  loop  occurs.  This  equals  one 
iteration  of  the  interesting  code. 

time_for_one_iteration  := 

(  my_elapsed_time  -  dummy_elapsed_time  )  /  loop_count  ); 

note  that  the  value  of  my_timer  is  not  changed  by  a  call  to 
timer. elapsed_tirae.  You  can  have  as  many  timers  as  you  wish. 
Revision  history;  S.  French  25-Jan-1989 
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package  body  timer  is 


—  Package:  timer 


—  Description:  This  package  provides  two  subprograms  which  are  called 

—  to  time  functions  in  the  unix  environment.  This  timer  was  used  to 

—  perform  timing  on  the  MIPS  Magnum  3000  computer  system. 


>  —  References:  Please  see  the  man  page  clock  for  information  about  how 

—  the  Unix  time  is  obtained. 


—  Revision  history:  S.  French  25-Jan-1989 


—  Host  processor:  gppe  gppe  3_10  UMIPS  mips  ml20-5  ATT  V3  0 


—  Package  dependencies:  This  package  interfaces  with  the  UNIX  system 

—  call  "clock". 


function  unix_clock  return  integer; 
pragma  INTERFACE(C,  unix_clock); 
pragma  INTERFACE_NAME(unix_clock, "clock" ) ; 

procedure  init_timer(  obj  :  out  raicro3ec_timer  )  is 
begin 

obj  :*  microsec_timer(  unix_clock  ); 
end  init_timer; 

function  elapsed_time(  obj  :  in  microsec_timer  )  return  Integer  is 
begin  ~ 

return  unix_clock  -  integer(  obj  ); 
end  elapsed_time; 

function  identity(  arg  :  integer  )  return  integer  is 
sorae_value  :  integer  :■=  0; 
begin 

some_value  :=  some_value  +  arg; 
return  some_value; 
end  identity; 

function  always_true  return  boolean  is 
bool_value  :  boolean  :=>  true; 
begin 

return  bool_value; 
end  always_true; 

begin 

null; 

end  timer; 
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with  text_io;  use  textio; 

with  data_struct_pkg;  use  data_struct_pkg; 

procedure  READDATA  (INFILE  :  IN  OUT  F1LE_TYPE; 

positions  :  out  means_array_type; 

weights  ;  out  weights_array_type)  is 


—  Procedure  Name:  READ_DATA 

—  Purpose:  Procedure  READ_DATA  is  called  by  Procedure  Main 

to  read  in  the  input  data  file. 


i  :  integer  ;=  0; 
j  :  integer  :«  0; 

package  FL0AT_I0  is  new  TEXT_I0. FL0AT_I0( FLOAT) ; 
begin 

while  not  END_OF_FILE(INFILE)  loop 
FL0AT_I0 . GET ( INFI LE ,  pos i t i ons ( i ) ) ; 
i  :»  i  +  1; 

FL0AT_I0 . GET< INFILE ,  pos i t ions ( i )) ; 
i  :=  i  +  1; 

FL0AT_I0.GET( INFILE,  weights(j)); 
j  j  1; 

SKIP_LINE( INFILE) ; 
end  loopj 

CLOSE(INFILE); 

end  READ  DATA; 
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58.000000  19.000000  62.000000 
59.000000  19.000000  74.000000 
60.000000  19.000000  70.000000 
57.000000  20.000000  76.000000 
58.000000  20,000000  90.000000 
59.000000  20.000000  91.000000 
60.000000  20.000000  78.000000 
61.000000  20.000000  52.000000 
56.000000  21.000000  75.000000 
57.000000  21.000000  96.000000 
58.000000  21.000000  101,000000 
59.000000  21.000000  91.000000 
60.000000  21.000000  67.000000 
56.000000  22.000000  75.000000 
57.000000  22.000000  100.000000 
58.000000  22.000000  108.000000 
59.000000  22.000000  104.000000 
60.000000  22.000000  92.000000 
61.000000  22.000000  72.000000 
56.000000  23.000000  84.000000 
57.000000  23.000000  135.000000 
58.000000  23.000000  159.000000 
59.000000  23.000000  165.000000 
60.000000  23.000000  155.000000 
61.000000  23.000000  124.000000 
62.000000  23.000000  61.000000 
56.000000  24.000000  105.000000 
57.000000  24,000000  164.000000 
58.000000  24.000000  190.000000 
59.000000  24.000000  196.000000 
60.000000  24.000000  185.000000 
61.000000  24.000000  151.000000 
62.000000  24.000000  74.000000 
55.000000  25.000000  59.000000 
56.000000  25.000000  109,000000 
57.000000  25.000000  173.000000 
58.000000  25.000000  200.000000 
59.000000  25.000000  206.000000 
60.000000  25.000000  194.000000 
61.000000  25.000000  160.000000 
62.000000  25.000000  81.000000 
63.000000  25.000000  58.000000 
36.000000  26.000000  52.000000 
56.000000  26.000000  105.000000 
57.000000  26.000000  165.000000 
58.000000  26.000000  191.000000 
59.000000  26.000000  197.000000 
60.000000  26.000000  185.000000 
61.000000  26,000000  152.000000 
62.000000  26.000000  80.000000 
63.000000  26.000000  58.000000 
36.000000  27.000000  60.000000 
37.000000  27.000000  73.000000 
40.000000  27.000000  85.000000 
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